A  PROGRAMMABLE  LIQUID  COLLIMATOR  FOR  BOTH  CODED 
APERTURE  ADAPTIVE  IMAGING  AND  MULTIPLEXED  COMPTON 

SCATTER  TOMOGRAPHY 

THESIS 

Jack  G.  M.  FitzGerald,  2d  Lt,  USAF 
AFIT/NUCL/ENP/12-M01 

DEPARTMENT  OF  THE  AIR  FORCE 
AIR  UNIVERSITY 

AIR  FORCE  INSTITUTE  OF  TECHNOLOGY 

Wright- Patterson  Air  Force  Base,  Ohio 


DISTRIBUTION  STATEMENT  A 

APPROVED  FOR  PUBLIC  RELEASE;  DISTRIBUTION  UNLIMITED 


The  views  expressed  in  this  document  are  those  of  the  author  and  do  not  reflect  the 
official  policy  or  position  of  the  United  States  Air  Force,  the  United  States  Department 
of  Defense  or  the  United  States  Government.  This  material  is  declared  a  work  of  the 
U.S.  Government  and  is  not  subject  to  copyright  protection  in  the  United  States. 


AFIT/NUCL/ENP/12-M01 


A  PROGRAMMABLE  LIQUID  COLLIMATOR  FOR  BOTH  CODED 
APERTURE  ADAPTIVE  IMAGING  AND  MULTIPLEXED  COMPTON 

SCATTER  TOMOGRAPHY 


THESIS 


Presented  to  the  Faculty 
Department  of  Engineering  Physics 
Graduate  School  of  Engineering  and  Management 
Air  Force  Institute  of  Technology 
Air  University 

Air  Education  and  Training  Command 
in  Partial  Fulfillment  of  the  Requirements  for  the 
Degree  of  Master  of  Science 


Jack  G.  M.  FitzGerald,  BS 
2d  Lt,  USAF 


March  2012 


DISTRIBUTION  STATEMENT  A 

APPROVED  FOR  PUBLIC  RELEASE;  DISTRIBUTION  UNLIMITED 


AF  IT/M ’CL  ENP.  12-M01 


A  PROGRAMMABLE  LIQUID  COLLIMATOR  FOR  BOTH  CODED 
APERTURE  ADAPTIVE  IMAGING  AND  MULTIPLEXED  COMPTON 

SCATTER  TOMOGRAPHY 


lack  G.  M.  FitzGerald,  BS 
2d  Lt.  USA  F 


Approved: 


2-14/  ?2*4LA 

Larry  W/feurggraf,  PhD  ( 


Benjamin  R  Kowash  (Member) 


Date 


I  £ 

Date 

Date 


AFIT/NUCL/ENP/12-M01 


Abstract 

A  novel,  fully  reconfigurable  collimator  device  for  7-ray  and  X-ray  imaging  was  built 
and  tested  as  a  coded  aperture.  The  device  consisted  of  10x10,  5x5x5  mm3  cham¬ 
bers.  Each  chamber  either  was  filled  with  an  attenuating  liquid,  stopping  photons, 
or  evacuated  of  the  attenuating  liquid,  allowing  the  photons  to  pass  through.  As  the 
pattern  of  “on”  and  “off”  chambers  was  manipulated,  different,  semi-independent 
views  of  the  7-ray  source  were  found.  Noise  in  reconstructed  images  decreased  in  all 
tests.  Image  reconstruction  was  performed  with  correlation  methods  and  Maximum 
Likelihood  Expectation  Maximization  (ML-EM).  With  ten  mask  patterns,  the  signal- 
to-noise  ratio  (SNR)  in  images  of  a  Co-57  point  source  increased  by  a  factor  of  4.3 
using  correlation  methods  and  by  a  factor  of  at  least  50  using  ML-EM.  SNR  in  images 
of  a  Cd-109  source  with  high  background  increased  by  a  factor  of  3.0  using  correlation 
methods  and  by  a  factor  of  1.8  with  ML-EM.  Two  extended  sources  were  imaged,  and 
the  images  improved  when  more  masks  were  used.  The  Multiplexed  Compton  Scat¬ 
ter  Tomography  (MCST)  forward  problem  using  a  PHDs  Co  high  purity  germanium 
(HPGe)  detector  was  tested  and  evaluated.  Potential  applications  are  discussed  in 
detail. 
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A  PROGRAMMABLE  LIQUID  COLLIMATOR  FOR  BOTH  CODED 
APERTURE  ADAPTIVE  IMAGING  AND  MULTIPLEXED  COMPTON 


SCATTER  TOMOGRAPHY 

I.  Introduction 


1.1  Overture 

Numerous  scientific  and  engineering  endeavors  rely  on  the  ability  to  image  ra¬ 
dioactive  sources  that  emit  high-energy  photons,  either  X-rays  or  7-rays.  Detection 
of  fissile  material,  nondestructive  material  testing,  medical  imaging  and  astronomy 
are  a  few  notable  examples.  This  paper  focuses  on  7-rays,  although  all  of  the  prin¬ 
ciples  and  techniques  are  equally  valid  for  X-rays,  particularly  hard  X-rays.  7-ray 
imaging  is  difficult  for  two  main  reasons.  First,  significant  depth  is  required  to  stop  a 
given  photon  fully  in  order  to  measure  the  energy  deposited.  The  mean  free  path  in  a 
detector  can  be  longer  than  the  detector’s  depth,  resulting  in  low  detection  efficiency. 
The  charge  displaced  by  the  incoming  photons  must  be  collected  by  equipment  that 
can  be  bulky,  increasing  the  overall  size  of  the  device.  7-ray  detectors  never  will  have 
the  pixel  density  of  a  visible  light  detector.  Secondly,  high-energy  photons  cannot 
effectively  be  diffracted  and  focused  by  a  lens.  The  wavelength  is  many  orders  of 
magnitude  smaller  than  the  characteristic  sizes  of  a  detection  device.  For  photons 
above  about  30  keV,  diffraction  may  be  ignored  [27,  45]. 

Since  the  photons  cannot  be  focused,  they  must  be  constrained  in  some  other  way. 
Often  a  pinhole  collimator  is  used.  A  pinhole  collimator  is  simply  a  small  hole  drilled 
into  a  sheet  of  high-Z  material.  When  the  pinhole  collimator  is  placed  between  a  7-ray 
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source  and  a  position-sensitive  detector,  all  of  the  7-rays  are  blocked  except  for  those 
that  pass  through  the  pinhole.  An  inverted  image  is  projected  onto  the  detector.  As 
the  pinhole  diameter  decreases,  the  spatial  resolution  of  the  system  increases.  The 
image  becomes  sharper.  However,  decreased  pinhole  diameter  also  causes  decreased 
throughput.  A  smaller  pinhole  corresponds  to  increased  detection  durations  since  a 
minimum  number  of  detected  photons  are  needed  to  form  a  good  image. 


* 


SOURCE 


RECONSTRUCTED 
IMAGE 


DECODING 
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Figure  1.  Left:  A  pinhole  collimator.  Right:  Basic  overview  of  the  coded  aperture 
imaging  procedure  [26]. 


One  way  to  increase  throughput,  while  maintaining  good  image  resolution,  is  to 
drill  more  holes  in  the  sheet  of  high-Z  material.  This  is  called  a  coded  aperture.  The 
raw  data  from  the  detector  consists  of  the  projections  from  each  of  the  pinholes,  all 
summed  together.  Knowing  the  positions  of  the  pinholes  relative  to  the  detector, 
a  computer  routine  can  be  used  to  unfold  all  of  the  projections  and  reconstruct  an 
image  of  the  source  [21,  2],  The  coded-aperture  technique  has  been  used  successfully 
for  many  decades.  The  main  drawback  to  the  technique  is  that  no  computer  routine 
can  reconstruct  an  image  of  the  source  perfectly.  All  reconstruction  routines  introduce 
what  is  termed  inherent  noise  to  the  final  image.  Some  inherent  noise  is  related  to 
the  pinhole  positions  whereas  some  inherent  noise  is  randomly  placed.  In  either  case, 
if  different  mask  patterns  are  used,  and  the  reconstructed  images  from  the  different 
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patterns  are  summed  together,  the  inherent  noise  will  cancel  out.  Therein  lies  the 
purpose  of  the  programmable  collimator.  By  reconfiguring  the  coded-aperture  during 
data  collection,  the  inherent  noise  cancels  out  in  the  final  image. 

In  2010,  at  the  Air  Force  Institute  of  Technology  (AFIT),  Dr.  Larry  Burggraf 
conceived  of  the  programmable  aperture  collimator  concept.  Confirmation  of  a  pro¬ 
grammable  collimator’s  feasibility  and  utility  was  provided  by  Maj  Benjamin  Kowash. 
They  pictured  a  device  that  would  use  a  liquid,  high-Z  metal  to  block  7-rays,  jux¬ 
taposed  with  a  low-Z  plastic  through  which  7-rays  could  pass.  If  the  high-Z  metal 
and  the  low-Z  plastic  could  be  moved  around  relative  to  each  other,  incident  7-rays 
would  be  blocked  in  different  patterns.  To  achieve  this  goal,  a  grid  was  drilled  out  of  a 
5x5x5  cm3  steel  block.  Small  plastic  plugs  were  fabricated  to  snugly  fit  in  to  each  grid 
chamber.  The  plugs  are  pushed  or  pulled  back  and  forth  in  the  chambers.  From  one 
side  of  the  grid  block,  the  liquid,  high-Z  metal  is  allowed  to  flow  freely.  By  pushing 
the  plug  backward,  the  high-Z  liquid  metal  is  displaced  from  the  chamber,  leaving 
only  air  and  plastic  to  attenuate  a  given  7-ray.  If  the  plug  is  pulled  forward  in  a  given 
chamber,  the  chamber  fills  with  the  liquid  metal  so  that  7-rays  are  blocked  in  that 
particular  chamber.  In  effect  this  setup  can  be  thought  of  as  a  10x10  matrix,  where 
each  element  of  the  matrix  can  be  turned  either  “on”  or  “off”  to  7-ray  transmission. 

Two  potential  uses  for  this  device  are  presented  in  this  document: 

•  Coded- Aperture  Imaging:  As  mentioned  above,  a  coded-aperture  provides  good 
throughput  and  resolution  comparable  to  a  pinhole  collimator.  The  computer 
routine  introduces  noise  into  the  reconstructed  image.  By  using  different  pinhole 
configurations,  or  mask  patterns,  the  noise  and  artifacts  from  the  reconstruction 
process  cancel  out  in  the  final  image. 

•  Multiplexed  Compton  Scatter  Tomography:  A  radionuclide  and  a  7-ray  detec¬ 
tor  are  placed  on  one  side  of  a  sample.  7-rays  are  Compton  scattered  within 


3 


a  material  of  interest.  The  scattered  7-rays  are  collected  by  a  position  and  en¬ 
ergy  sensitive  detector.  An  image  of  the  electron  density  in  the  sample  is  found, 
which  is  useful  in  nondestructive  testing  or  medical  imaging,  especially  bone 
density  studies.  The  problem  is  ill-conditioned,  and  could  benefit  from  con¬ 
straints  through  physical  collimation.  Current  commercial  systems  are  highly 
collimated  and  examine  one  voxel  in  the  sample  at  a  time.  By  viewing  mul¬ 
tiple  voxels  simultaneously,  acquisition  time  decreases  and/or  required  source 
activity  decreases. 


1.2  Coded  Aperture  Application 

The  coded  aperture  system  was  first  suggested  by  Dicke  [21]  and  Abies  [2]  in  1968. 
X-ray  astronomers  had  been  seeking  a  system  in  which  heavenly  X-ray  sources  could 
be  imaged  with  good  resolution  and  in  the  shortest  time  possible.  A  pinhole  collimator 
is  used  to  increase  the  resolution  of  an  X-ray  camera  system.  The  smaller  the  pinhole 
diameter  is,  the  better  the  resolution  will  be.  However,  with  a  smaller  pinhole  the 
number  of  counts  incident  on  the  detector  decreases.  Dicke  and  Abies  suggested 
that  an  array  of  pinhole  collimators  could  be  placed  in  front  of  the  detector.  The 
pinholes  are  small  so  as  to  provide  good  resolution.  The  greater  quantity  of  pinholes 
increases  the  count  rate.  This  pinhole  array,  or  coded-aperture,  can  be  configured 
in  a  variety  of  ways.  As  light  from  a  source  travels  through  the  pinholes,  it  is  cast 
onto  the  detector  in  a  different  way  from  each  pinhole.  Each  pinhole  causes  a  unique 
projection  onto  the  detector.  A  computer  must  be  used  to  interpret  the  data  from  the 
detector,  unfolding  all  the  pinhole  projections  to  arrive  at  an  image  that  reasonably 
recreates  the  original  object.  A  coded  aperture  increases  the  Signal-to-Noise  Ratio 
(SNR)  of  a  standard  detector  [21],  The  SNR  advantage  of  the  coded- aperture  scheme 
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is  most  prominent  with  point  sources  and  becomes  less  pronounced  as  more  sources 
are  placed  in  the  Ficld-of-View  (FOV). 

Multiple  coded-aperture  mask  configurations  have  been  used  to  increase  the  fi¬ 
delity  of  a  reconstructed  image  [16].  Alternatively,  the  mask  and/or  detector  can 
be  moved  to  different  locations  around  the  object  to  increase  the  number  of  unique 
views.  Prefabricated  masks  have  been  used  for  7-ray  coded-apertures,  hitherto.  Of¬ 
ten  tungsten  is  used,  or  other  heavy  metals  such  as  lead.  For  7-rays  in  the  100’s  to 
1000’s  of  keV,  a  coded-aperture  thickness  on  the  order  of  centimeters  is  required.  The 
closed  portions  of  the  coded-aperture  collimator  must  be  thick  enough  to  attenuate 
a  large  portion  of  the  7-rays.  The  programmable  collimator  of  this  study  is  a  10x10 
array,  where  each  array  element  can  be  turned  “on,”  meaning  that  7-rays  are  allowed 
through,  or  “off,”  meaning  that  7-rays  are  attenuated.  The  ability  to  configure  the 
collimator  in  different  ways  allows  for  multiple,  semi-independent  views  to  be  made 
of  the  object,  increasing  image  fidelity  over  a  given,  total  measurement  time. 

As  a  dynamic  coded-aperture,  the  programmable  collimator  has  many  potential 
applications  in  the  Department  of  Defense  (DoD)  and  elsewhere.  These  applications 
are  discussed  further  in  Sections  2.6.5  and  6.2.6. 

•  The  programmable  collimator  could  be  mated  with  a  position  and  energy  sen¬ 
sitive  detector  and  employed  as  an  imaging  7-ray  detector  on  an  unmanned, 
ground-based  robot  or  aerial  vehicle.  The  SNR  advantage  from  the  coded- 
aperture  technique  would  be  advantageous  in  environments  where  a  source  is 
either  weak  or  shielded.  A  vehicle  of  this  type  could  survey  large  areas  after  a 
nuclear  attack  or  a  nuclear  reactor  accident,  as  well  as  find  certain  radioactive 
nuclides  in  effluents  for  treaty  monitoring  applications. 

•  A  passive  sensor  could  continuously  search  for  radioactive  material  in  passing 
vehicles  at  ports,  gates,  etc. 
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•  The  lessons  learned  from  programmable  collimation  could  be  applied  to  surveil¬ 
lance  satellites.  A  future  system  might  consist  of  a  number  of  detectors,  placed 
on  one  satellite  or  several  satellites,  watching  the  Earth  for  nuclear  events. 
The  same  research  conducted  for  choosing  the  best  mask  sequence  for  the  pro¬ 
grammable  collimator  is  applicable  to  space  surveillance  systems. 

•  The  device  could  allow  for  more  effective  7-ray  or  X-ray,  space-based  astron¬ 
omy.  Coded-apertures  have  a  long  history  with  astronomy.  Sometimes  SNR 
is  increased  by  spinning  the  satellite.  Programmable  collimation  removes  the 
need  to  spin,  since  spinning  can  introduce  difficulties  in  satellite  design. 

•  Medical  tomographical  imaging  could  be  improved  by  increasing  image  quality 
and/or  decreasing  the  radioactive  dose  given  to  patients. 

•  The  device  could  be  used  to  give  position  sensitivity  to  a  position-insensitive 
detector.  Such  functionality  has  been  explored  with  Rotating  Mask  Collimators 
(RMC)  mated  with  Sodium  Iodide  (Nal)  detectors  [74],  A  programmable  col- 
limator/Nal  system  would  be  cheaper  than  most  energy  and  position  sensitive 
detectors.  Such  a  system  could  be  used  in  any  of  the  applications  mentioned 
above,  among  others. 

1.3  Multiplexed  Compton  Scatter  Tomography  Application 

Another  type  of  application  for  the  programmable  collimator  is  the  improvement 
of  Multiplexed  Compton  Scatter  Tomography  (MCST).  MCST  was  first  investigated 
at  AFIT  by  Brian  Evans  and  Matt  Lange  in  research  advised  by  Professors  Jeffrey 
Martin  and  Larry  Burggraf  [25,  41],  That  research  demonstrated  multiplexed  imaging 
in  a  single  discretized  slice  of  material.  MCST  is  a  methodology  for  one-sided  non¬ 
destructive  investigation  of  materials.  It  differs  from  other  tomographic  techniques 
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because  the  source  and  detector  can  be  placed  on  the  same  side  of  the  object,  rather 
than  on  opposite  sides.  In  many  industrial  cases,  access  to  opposite  sides  of  the  object 
to  be  imaged  is  not  available.  MOST  could  be  applied  effectively  toward  the  detection 
of  cracks  and  structural  failures  in  airfoils,  composite  structures,  and  other  aircraft 
systems.  MOST  also  could  be  used  for  improved  medical  tomography.  Probably 
it  would  complement  existing  tomographic  techniques,  like  Computed  Tomography 
(CT)  or  Positron  Emission  Tomography  (PET),  allowing  for  better  resolution  and/or 
a  lower  total  dose  to  patients.  In  MOST,  a  monoenergetic  radionuclide  is  placed  near 
an  object,  beside  a  position  and  energy  sensitive  detector.  Shielding  is  placed  between 
the  source  and  detector.  Photons  emitted  by  the  source  are  Compton  scattered  in 
the  sample.  Some  of  these  scattered  photons  are  detected.  By  knowing  the  position 
and  energy  of  the  received  photons,  spatial  and  compositional  information  about  the 
sample  is  ascertained.  MCST  is  well  suited  for  two-dimensional  applications.  The 
programmable  collimator  would  serve  to  produce  three-dimensional  images  by  con¬ 
straining  the  highly  multiplexed  problem,  allowing  for  solutions  to  be  found  quickly 
and  accurately. 

1.4  Research  Goals 

The  goal  of  this  research  was  the  characterization  of  the  programmable  collimator 
given  the  coded-aperture  and  MCST  setups.  To  the  best  knowledge  of  the  author 
and  his  committee,  a  liquid  metal  never  has  been  used  to  attenuate  7-rays  in  quickly 
programmable  patterns.  An  evaluation  of  the  viability  of  this  device  for  different 
applications  was  wanted.  Recommendations  from  this  proof-of-concept  study  will  be 
used  to  direct  future  investigations  into  programmable  collimation. 
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II.  Theory  and  Background 


2.1  Coded- Aperture  Theory 

The  coded  aperture  system  was  first  suggested  by  Dicke  and  Abies  in  1968  [21,  2], 
In  7-ray  imaging,  a  single  pinhole  collimator  is  used  to  increase  the  resolution  of  a 
detector.  With  coded  apertures  the  single  pinhole  is  replaced  by  numerous  pinholes 
in  an  array.  As  high-energy  light  shines  upon  the  array,  each  pinhole  causes  a  unique 
projection  onto  a  spatially-sensitive  detector.  Knowing  the  configuration  of  the  pin¬ 
hole  array,  mathematical  techniques  are  used  to  unfold  all  the  projected  images  from 
each  other,  reconstructing  the  original  image.  The  term  multiplex  advantage  is  used 
to  describe  the  greater  SNR  in  reconstructed  images  produced  by  a  coded-aperture 
array  versus  those  from  a  single  pinhole  (SNR  is  explained  thoroughly  in  Section 
5.3).  For  point  sources  the  SNR  improvement  goes  roughly  as  \/N,  where  N  is  the 
number  of  pinholes  [21].  As  more  point  sources  are  imaged,  the  multiplex  advantage 
decreases.  A  coded-aperture  works  better  for  point  sources  than  for  extended  sources. 

Coded- aperture  problems  usually  are  approximated  as  three  parallel  planes.  Fen- 
imore  and  Cannon  outlined  a  basic  coded-aperture  methodology  in  their  1978  paper 
[26].  The  sources  exist  in  the  object  plane,  O.  The  radiation  moves  through  a  col¬ 
limator  in  the  aperture  plane,  A,  where  A  is  the  mask  pattern.  The  light  moving 
through  the  aperture  plane  is  projected  onto  the  detector  or  picture  plane,  P,  given 
by 

P=(0*A)  +  N  (1) 

As  with  all  real  detector  experiments,  noise,  N,  always  is  summed  with  the  true 
signal  and  must  be  considered.  Sources  of  noise  include,  but  are  not  limited  to,  low 
counting  statistics,  background  scattering,  transmission  through  closed  portions  of 
the  array,  imprecisions  in  the  detector  and  electronics  noise.  In  most  applications,  all 


three  planes  are  discretized  into  two-dimensional  matrices.  Notice  that  O  is  convolved 
with  A  in  Equation  1,  where  *  is  the  convolution  operator.  An  estimate  of  the  object, 
O,  can  easily  be  found  as 
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where  J-  is  the  Fourier  transform  and  R  is  the  reflection  operator.  This  method  is 
typically  flawed,  because  T  (A)  contains  regions  of  frequency  with  small  magnitudes. 
In  these  regions  the  true  signal  easily  is  lost  in  the  noise.  Woods  et  al.  showed 
that  deconvolution  techniques  can  be  used  in  conjunction  with  a  two-dimensional 
Wiener  filter  [77].  Often  correlation  techniques  are  preferred.  Here  *  is  the  correlation 
operator.  Matched  filtering  is  another  term  used  to  describe  correlation  decoding. 

In  the  correlation  methodology,  the  object  reconstruction  is  defined  as 


6  =  P*G 

=  RO  *  (A*G)  +  N  *G 


(3) 


G  is  called  the  postprocessing  array,  or  matched  filter,  and  is  only  a  mathematical 
construct.  It  does  not  exist  physically.  G  is  chosen  such  that  A  *  G  approximates  a 
two-dimensional  Dirac  5-function.  If  this  5-function  were  perfect,  the  object  estimate 
would  reduce  to 

0  =  0  +  (N*G)  (4) 

Notice  that  the  noise  cannot  be  removed  from  the  estimated  object.  Unfortunately, 
A-kG  will  deviate  from  a  true  5-function  in  any  real  application,  introducing  artifacts 
into  the  estimated  object. 

In  the  original  astronomy  application,  a  point  source  is  assumed  to  be  infinitely 
far  away.  A  star  is  much  more  distant  than  the  scale  of  the  detector  apparatus.  If  the 
source  is  at  infinity,  the  aperture  pattern  is  projected  onto  the  detector  plane  at  the 
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Figure  2.  Planar  approximation  of  the  coded-aperture  problem.  Adapted  from  [65]. 


same  scale.  However,  in  most  defense  and  medical  applications,  the  source  cannot  be 
assumed  to  be  located  at  infinity.  Rather,  a  magnification  effect  will  exist,  and  must 
be  considered.  The  magnification  factor,  m,  is  given  as: 

-Sl  +  S2 

m  =  - 

-si 

where  .s  i  is  the  distance  between  the  object  and  the  aperture  planes  and  S2  is  the 
distance  between  the  aperture  and  detector  planes,  as  shown  in  Figure  2.  In  near-field 
applications,  it  is  conceivable  that  the  magnification  could  be  used  to  determine  the 
distance  of  the  source  from  the  detector.  Information  about  the  source’s  location  in  all 
three  dimensions  could  possibly  be  gathered  by  factoring  in  the  effect  of  magnification. 


2.2  Mask  Configurations 

A  number  of  different  coded-aperture  collimator  designs  have  been  employed  in  the 
past,  each  with  its  own  advantages.  Fresnel  zone  plates,  annuli,  and  Fourier  apertures 
can  be  used  as  coded- apertures  [9].  The  programmable  collimator  falls  under  the 
broad  pinhole-array  category.  The  term  “pinhole”  is  used  here  somewhat  loosely  as 
any  small  round  or  square  aperture.  Each  element  of  the  programmable  collimator  is 
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shaped  as  a  5x5  mm2  hole,  for  instance.  Square  holes  are  easy  to  manufacture,  and 
thus  are  fairly  common.  Hexagonal  holes  also  have  been  organized  into  a  honeycomb 
pattern  and  used  successfully  [17]. 

The  first  pinhole  array  design  considered  is  called  the  random  array  and  was 
suggested  in  Dicke’s  original  paper  [21].  Any  number  of  pinholes  can  be  opened,  with 
varying  results.  Often  50%  of  the  array  elements  are  opened,  as  is  the  case  in  Figure 
3.  An  initial  guess  for  the  matched  filter,  G,  could  be  the  matrix  A  itself.  Since  a 
source  projects  an  identical  shadow  of  the  aperture  pattern  onto  the  detector,  this 
is  a  reasonable  first  guess.  However,  recall  that  the  term  A  *  G  should  approach  a 
5-function  as  nearly  as  is  possible.  If  A  is  used  for  G,  a  distinct  triangular  base  is 
formed  in  the  correlation  (Figure  3).  Anytime  that  the  sidelobes  of  the  correlation 
deviate  from  zero,  artifacts  are  introduced  into  the  object  estimate,  O. 


Autocorrelation  of  the  random  pinhole  array,  A 


Slice  of  the  autocorrelation  of  A. 


Figure  3.  Top:  A  random  pinhole  array,  A.  Bottom-Left:  The  autocorrelation  of  the 
pinhole  array,  A.  Bottom-Right:  A  slice  in  one  dimension  of  the  autocorrelation  of  A. 
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Of  course,  G  is  nothing  but  a  mathematical  abstraction,  and  can  be  altered  in 
any  way  desired.  With  this  idea  in  mind,  C.  Brown  suggested  that  the  anti-mask 
should  be  allowed  to  take  on  negative  values  [10].  Namely,  the  zeros  in  the  mask  are 
replaced  with  “-l’s”,  or 


G{i,j)  =  1  if  A(i,j)  —  1 
=  — 1  if  A(i,j)  —  0 


(6) 


where  i  and  j  are  row  and  column  indices,  respectively.  In  the  mask,  A,  a  value  of  “1” 
corresponds  to  an  open  element  where  y-rays  are  allowed  through,  whereas  a  value 
of  “0”  corresponds  to  closed  elements,  where  y-rays  are  blocked. 


Correlation  of  the  mask.  A,  with  a  Brown  type  anti-mask,  G 


Slice  of  the  correlation  of  the  random  array  mask,  A,  with  a  Brown  type  anti-mask,  G 


Figure  4.  Left:  The  correlation  of  A  from  Figure  3,  with  its  Brown  anti-mask,  G.  G  is 
constructed  by  replacing  the  zeros  in  A  with  “-l’s”.  Right:  A  2D  slice  of  the  3D  AkkG 
plot. 


The  Brown  anti-mask  improves  the  A  *  G  term.  The  triangular  base  is  removed, 
causing  the  term  to  approximate  a  5-function  more  closely.  Note  that  the  inclusion  of 
negative  values  in  the  anti-mask  causes  negative  values  in  the  A  *  G  array.  This  can 
result  in  non-physical  artifacts  in  O,  though  the  advantage  of  the  removed  triangular 
base  typically  outweighs  the  disadvantage  of  these  artifacts. 
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Another  mask  configuration  is  called  the  Uniformly  Redundant  Array  (UR A). 
A  subset  of  this  mask  class  is  the  Modified  Uniformly  Redundant  Array  (MURA). 
URA’s  were  studied  extensively  by  Fenimore  and  Cannon  [26].  In  an  extension  of 
Calabro  and  Wolf’s  work,  they  define  the  aperture  plane  as  an  infinite  mosaic  of 
repeating,  pseudonoise  arrays  [12].  The  basic  array  has  dimensions  of  r  and  s  where 
r  and  s  are  both  prime  numbers  and  r  —  s  +  2.  The  basic  array,  A(I,  J)  is  defined  as 


A(I,  J )  =  0  if  /  =  0, 

=  1  if  J  =  0  and  7^0, 
=  1  if  Cr{I)Cs(J)  =  1, 
=  0  otherwise 


(7) 


where  7  and  J  are  row  and  column  indices  of  the  basic  array,  respectively,  and 


Ca(f3)  =  1  if  there  exists  an  integer  x,  1  <  x  <  a 

such  that  7  =  raodQr2  (8) 

=  —1  otherwise 


where  a  can  be  r  or  s  and  f3  can  be  7  or  J.  Physical  coded-apertures  cannot  be 
infinite,  but  instead  are  a  finite  cutout  of  the  infinite  mosaic.  Figure  5  shows  a  coded- 
aperture,  A,  of  size  2 r  x  2s,  where  r  =  19  and  s  =  17.  G,  the  Brown  matched  filter, 
was  found  as  in  Equation  6. 

A  discrete  correlation  always  doubles  the  size  of  the  output  matrix  as  compared 
to  the  input  matrices.  The  center  square  of  the  correlation  is  the  only  part  used  in 
the  image  reconstruction  process.  The  additional  spikes  in  (A  *  G )  arise  from  the 
repetition  of  the  basic  array  within  A.  The  URA  has  been  shown  to  outperform  the 
random  array  in  most  cases,  when  correlation  methods  are  used  [5,  26].  Unfortunately, 
the  programmable  collimator  is  not  equipped  with  sufficient  elements  for  a  proper 
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URA/MURA  pattern.  Technically  one  conld  be  formed,  but  whether  the  resulting 
matrix  could  be  classified  as  a  URA  would  be  disputable.  As  such,  random  arrays  will 
be  focused  upon  with  the  programmable  collimator.  On  top  of  all  this,  URA/MURAs 
are  not  necessarily  more  effective  for  ML-EM  methods. 


Example  Unifoimly  Redundant  Airay  MS,  5=17 
Wlite- Open  toy's  Blacks  Closed  toy's 


y 


Figure  5.  Top:  An  example  URA,  A,  where  r=19,  s=17  and  A  is  2r  x  2s.  Bottom-Left: 
The  correlation  of  the  URA  mask,  A,  with  the  Brown  anti-mask,  G.  Bottom-Right: 
The  central  slice  of  the  correlation  between  A  and  G.  [38] 

Note  that  the  above  figures  do  not  have  any  noise  included.  The  portions  in  the 
plots  of  (A-kG)  that  appear  to  be  noise  are  instead  artifacts  of  the  correlation  process. 
If  the  A  and  G  planes  approach  infinite  extent  in  space,  the  artifacts  would  diminish 
to  zero.  These  artifacts  often  are  termed  inherent  noise.  The  choice  of  aperture 
pattern  and  matched  filter  pattern  should  be  made  so  as  to  diminish  the  inherent 
noise  as  much  as  possible. 


14 


2.3  Reconfiguring  the  Mask 


The  primary  advantage  of  the  programmable  collimator  is  its  ability  to  be  recon¬ 
figured.  Multiple  “views”  are  desired  of  the  imaged  object.  A  number  of  different, 
unique  projections  of  an  object  will  increase  the  fidelity  of  the  final,  reconstructed 
image.  This  can  be  achieved  in  a  variety  of  different  ways.  In  some  cases,  the  coded- 
aperture  mask  and/or  detector  is  rotated,  giving  multiple  views  of  the  object  [17]. 
Often,  in  this  case,  hexagonal  arrays  are  used,  being  more  conducive  to  rotational 
transformations  about  the  center  axis  of  the  mask.  In  the  satellite  application,  this 
sort  of  rotation  adds  complexity  to  the  satellite  system.  The  entire  satellite  must 
rotate,  perhaps  causing  problems  with  other  systems.  In  ground-based  applications, 
the  requirement  of  rotation  implies  additional  actuators  to  make  the  whole  system 
rotate. 

In  other  tests,  the  coded-aperture  mask  pattern  was  altered  by  sliding  a  pattern 
cut  from  tungsten  in  one  dimension  and  only  allowing  a  portion  of  the  pattern  to 
let  radiation  through  [71,  16].  For  example,  a  10  x  20  pattern  could  be  cut  from  the 
tungsten.  Lead  could  be  used  to  ensure  that  only  a  square  of  size  10  x  10  would  be 
open.  The  pattern  could  be  slid  in  the  long  direction,  allowing  for  ten  mask  patterns. 
Such  a  design  was  suggested  as  a  way  to  combat  the  incomplete  attenuation  of  photons 
in  the  closed  portions  of  a  coded-aperture  pattern,  which  becomes  a  greater  factor 
with  higher  energy  y-rays.  The  programmable  collimator  is  an  improvement  over 
this  sort  of  design  in  that  it  can  produce  any  possible  pattern,  therefore  allowing  for 
“views”  that  are  more  independent  from  each  other. 

Another  means  for  gaining  multiple  views  of  the  object  is  either  to  rotate  the 
object  or  to  rotate  the  coded- aperture/detector  system  around  the  object.  The  latter 
method  is  recognizable  as  a  staple  CT  technique.  This  sort  of  rotation  is  especially 
well  suited  for  three-dimensional  imaging,  or  tomography.  The  disadvantages  are 
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Figure  6.  A  time-coded  pinhole  array,  in  the  style  of  Clinthorne  et  al.  [16].  The 
columns  are  shifted  once  leftward  in  each  subsequent  array. 


similar  to  those  enumerated  previously,  namely  the  additional  mechanical  complexity 
of  a  rotating  system.  Furthermore,  the  rotation  can  cause  problems  for  iterative 
reconstruction  algorithms,  like  ML- EM.  Such  techniques  require  the  pixel  size  in  the 
object  be  fixed  so  that  the  relationship  between  the  object  and  the  projections  is 
tractable  [33]. 

With  the  programmable  collimator,  no  rotation  is  necessary.  Semi-independent 
views  are  gained  by  changing  the  mask  configuration.  Thus,  the  collimator  and 
detector  can  stay  in  the  same  position  relative  to  the  object.  Or  the  programmable 
collimator  could  be  conjoined  with  current  methods,  making  them  better. 

Though  not  investigated  in  this  work,  a  future  programmable  collimator  imaging 
system  might  be  able  to  incorporate  incoming  data  in  order  to  choose  the  best  patterns 
for  the  remainder  of  the  data  run.  Aperture  patterns  would  be  chosen  on  the  fly.  After 
the  first  1-2  arbitrary  mask  patterns  were  used,  the  system  would  have  a  good  idea  of 
where  the  source  is  located.  It  could  then  choose  the  proper  patterns  so  as  to  narrow  in 
on  the  source.  The  system  could  start  with  a  high  pinhole  density,  allowing  maximum 
throughput  with  lower  resolution.  As  the  source  becomes  visible,  the  pinhole  density 
would  decrease,  resulting  in  a  sharper  image  in  a  shorter  duration. 
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2.4  The  Maximum  Likelihood  Expectation  Maximization  Algorithm 

Matched  filter  and  correlation  techniques  are  the  traditional  means  for  analyzing 
data  from  a  coded-aperture  system.  A  more  modern  and  effective  method,  called  the 
Maximum  Likelihood  Expectation  Maximization  (ML-EM)  technique,  has  been  used 
for  several  decades  in  medical  imaging  [61].  In  ML-EM,  an  accurate  forward  model 
must  first  be  devised.  Good  predictions  of  the  detector  response  are  required,  given 
mask  position  and  configuration,  source  type,  source  location,  etc.  The  more  accurate 
the  forward  model  is,  the  better  the  results  given  by  the  ML-EM  method  will  be.  The 
forward  model  could  even  be  created  experimentally,  by  carefully  setting  up  a  source- 
mask-detector  system  and  measuring  all  of  the  relevant  parameters.  This  approach  is 
rarely  practical.  Instead,  computer  models  typically  are  used  as  the  forward  model, 
sometimes  incorporating  real-world  data. 

The  forward  model  is  used  to  define  the  ML-EM  variable  A.  The  forward  model 
deals  with  two  sets  of  spatial  locations,  being  the  detector  plane  and  the  source 
plane.  The  detector  is  discretized  into  N  pixels.  The  source  plane  is  discretized  into 
M  pixels.  In  more  advanced  versions,  the  source  region  could  be  split  into  three- 
dimensional  voxels.  In  this  study,  only  2D  pixels  are  considered.  The  purpose  of  the 
forward  model  is  to  predict  the  number  of  counts  that  would  be  recorded  by  a  given 
detector  pixel  if  the  source  were  located  at  a  given  source  plane  pixel.  The  model 
should  incorporate  as  many  physical  processes  as  possible,  such  as  inverse  squared 
geometric  attenuation,  material  attenuation  and  the  location  and  configuration  of 
the  coded-aperture  mask.  More  advanced  models  might  include  Compton  scatter, 
heterogeneous  efficiencies  among  detector  strips  (and  within  detector  strips),  and  so 
forth.  The  goal  of  the  forward  model  is  to  fill  in  an  NxM  matrix  named  A,  where  Aij 
is  the  predicted  response  of  detector  pixel  i  if  the  source  were  located  at  source  plane 
pixel  j  (See  Section  3.1).  The  values  of  the  elements  of  A  can  be  thought  of  loosely  as 


17 


the  expected  number  of  counts.  However,  A  must  be  normalized  before  being  used. 
The  most  appropriate  way  to  consider  Ai3  is  as  the  probability  that  a  photon  emitted 
from  a  source  at  location  j  is  detected  by  detector  pixel  i.  A  is  precomputed  and  is 
used  in  the  ML-EM  algorithm. 

Once  A  is  computed,  the  ML-EM  iterative  sequence  is  performed.  Besides  A,  the 
algorithm  requires  a  few  other  vectors,  y  is  a  length  N  column  vector  that  contains 
the  real  data  from  an  experiment.  That  is,  y%  is  the  number  of  counts  recorded  by 
detector  pixel  i.  A  is  the  object/source  estimate  and  is  the  output  of  the  algorithm.  A 
is  a  length  M  column  vector,  where  A  j  is  the  probability  that  the  source  was  located 
at  source  plane  pixel  j.  b  is  an  additive  noise  term  and  represents  the  probability  that 
a  detected  photon  is  attributable  to  background  or  Compton  scattering.  The  single 
equation  of  the  ML-EM  algorithm  is  [53,  40,  62]: 

\  old 

\new  _  j _  \  '  a  _ Vi _ 

J  i3EkAk\k  +  b 

In  MATLAB®  this  is  written  as: 


Anew  A 


Old 


A  *  <  y  ./  (A  *  A  + 


.ja 


(10) 


where  a  is  a  normalization  factor  for  A,  found  by  multiplying  the  inverse  of  A  by  a 
length  N  ones  vector.  To  perform  the  ML-EM  algorithm,  every  element  of  A  is  first 
set  to  “1.”  In  other  words  the  source  has  an  equal  probability  of  being  anywhere 
before  the  algorithm  begins.  The  code  then  continues  to  loop  Equation  9.  After  the 
equation  is  evaluated  in  each  loop,  A0id  is  set  to  Anew.  The  ML-EM  equation  compares 
the  real  data  to  the  predictions,  morphing  A  with  each  iteration  to  make  the  best 
agreement  between  the  prediction  and  reality.  There  is  no  set  way  to  determine  the 
minimum  number  of  iterations  necessary.  It  must  be  determined  empirically  for  any 
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given  problem.  It  is  possible  to  perform  too  many  iterations  (Section  5.2),  in  some 
cases.  Weak  regions  can  be  lost  when  too  many  iterations  are  performed.  The  user 
must  find  the  proper  number  of  iterations  in  order  to  reach  convergence  without  losing 
weak  regions. 

2.5  Multiplexed  Compton  Scatter  Tomography 

Multiplexed  Compton  Scatter  Tomography  (MOST)  has  been  explored  previously 
at  the  Air  Force  Institute  of  Technology  (AFIT).  Brian  Evans  and  Matt  Lange,  in 
research  advised  by  Professors  Jeffrey  Martin  and  Larry  Burggraf,  looked  at  defects 
in  thin  aluminum  airfoils  [25,  41].  Evans  used  six  HPGe  detectors  and  collimated  the 
scattered  photons  down  to  a  two-dimensional  slice.  MOST,  in  that  case,  would  be 
useful  in  detecting  cracks  that  develop  around  rivets  on  aircraft.  Marc  Sands  used 
MOST  to  image  a  phantom  representing  a  wrist  bone  [59].  Noninvasive  bone  density 
measurements  would  be  useful  for  the  medical  characterization  of  osteoporosis  in 
patients.  Sand’s  MOST  system  was  able  to  differentiate  between  normal,  osteoporotic, 
and  void  bone  densities. 

MOST  is  a  nondestructive  inspection  technique  whereby  an  object’s  near-surface 
interior  structure  may  be  discerned.  In  MOST,  high  energy  photons  shine  into  an 
object.  By  detecting  the  photons  that  Compton  scatter  back  toward  the  detector, 
information  about  the  elements  within  the  object  is  found.  This  technique  differs 
from  computed  tomography  (CT)  where  the  object  is  placed  between  the  source  and 
the  detector  and  the  attenuation  through  the  object  is  found.  With  MOST,  the  source 
and  the  detector  are  located  on  the  same  side  of  the  object.  The  CT-type  setup  is 
not  always  possible. 

MCST  takes  advantage  of  the  fact  that  the  energy  of  a  Compton  scattered  photon 
is  angle  dependent.  Furthermore,  the  cross  section  for  scatter  depends  on  the  electron 
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density,  and  thus  the  atomic  number,  of  the  scattering  element  in  question.  Compton 
scattering  is  governed  by  the  Klein-Nishina  equation  [37] 


da  7r2(  1  V  ( 1+^2N\  f  1  _| _ q2(!-/3)2  \ 

dfl  ^'0^1 +a(l-/3))  \  2  )  ^  "T"  (l+02)[l+a(l-0)] ) 

/3  =  cos  9 


(11) 


where  a  is  the  cross  section  for  scattering,  in  barns,  Tq  is  the  classical  electron  radius, 
9  is  the  scattering  angle,  D  is  the  solid  angle  into  a  given  detector  pixel,  and  a  = 
E/rrtQC2.  The  experimenter  is  most  interested  in  hireling  the  Z  of  the  material,  and 
thus  must  integrate  over  the  solid  angle,  which  is  often  9  dependent.  In  MOST,  the 
experimenter  gets  a  response  in  each  detector  pixel.  From  an  ideal  point  sample, 
this  response  would  be  a  ^-function  at  a  certain  energy,  multiplied  to  a  certain  count 
height.  The  energy  of  the  scattered  photon  is  given  by  the  Compton  formula  [37]: 


E' 


E 


1  + 


E 

rriQC2 


(1  —  cos  9) 


(12) 


where  E  and  E'  are  the  energy  of  a  photon  before  and  after  Compton  scattering, 
respectively.  By  looking  at  the  energy  of  the  detected  photon,  and  by  knowing  the 
initial  energy  from  a  nronoenergetic  source,  the  angle  of  scattering  can  be  deduced. 
The  number  of  counts  received  by  each  detector  pixel  would  also  be  known,  which  is 
directly  proportional  to  the  cross  section  for  scattering.  Given  all  this  knowledge,  the 
electron  density  of  the  material  point  is  found.  From  electron  density,  atomic  number 
can  follow. 

Note  that  the  Compton  equation  assumes  that  the  photon  scatters  off  of  an  ini¬ 
tially  stationary  electron.  In  reality  the  electrons  are  moving  about  the  atoms  at 
significant  speeds.  Depending  on  how  a  given  electron  is  moving  relative  to  the  de¬ 
tector,  the  scattered  energy  can  be  more  or  less  than  that  given  by  the  Compton 
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formula.  This  effect  is  called  Doppler  broadening,  and  serves  to  widen  the  detec¬ 
tor  response  peaks.  Doppler  broadening  in  MOST  problems  has  been  characterized 
previously  [25,  41]. 


Figure  7.  By  knowing  the  source  and  detector  positions,  as  well  as  the  initial  and 
scattered  energy  of  a  single  photon,  an  isogonic  arc  of  possible  scatter  locations  is 
found  [25]. 

From  a  measured,  Compton  scattered  photon,  an  infinitesimal  point  in  the  sample 
is  not  deduced,  but  rather  the  scattered  angle  into  the  detector.  For  an  individual, 
detected  photon,  only  an  isogonic  arc  of  possible  locations  in  the  sample  where  the 
scattering  angle  is  6  is  known.  Various  techniques  can  be  used  to  discern  the  point 
at  which  the  scattered  photon  originated,  rather  than  the  isogonic  arc.  Collimation 
can  limit  the  field  of  view  into  each  detector  pixel,  chopping  the  isogonic  arc  into 
small  pieces.  By  changing  the  collimation,  by  using  more  detectors,  and  by  changing 
the  source-object-detector  geometry,  the  true  point  of  scatter  can  be  determined 
from  all  the  different  isogonic  arcs.  The  isogonic  arcs  are  similar  to  the  straight  line 
projections  of  attenuation  in  computed  tomography.  Evans  and  Lange  observed  that 
the  electron  density  would  bleed  from  high  density  regions  to  low  density  regions  as  the 
reconstruction  routine  was  iterated.  The  line  between  high  and  low  electron  density 
regions  in  a  given  phantom  would  thus  appear  blurred  in  the  final,  reconstructed 
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image.  Collimation  causes  the  reconstruction  routine  to  be  more  stable,  reducing  this 
bleeding  effect. 

Perhaps  the  greatest  hindrance  to  MOST  is  the  fact  that  7-rays  are  attenuated  in 
the  object  before  and  after  being  scattered.  This  effect  alters  the  counted  events  in 
the  detector  as  compared  to  a  model  that  only  includes  Compton  scattering.  In  non¬ 
destructive  testing,  the  composition  of  an  object  being  imaged  is  obviously  not  known, 
meaning  that  the  attenuation  through  the  object  is  unknown.  Several  methods  have 
been  investigated  for  overcoming  the  problem  of  attenuation.  Lale  et  al.  ignored 
attenuation  altogether,  assuming  that  higher  energy  7-rays  (5.6  MeV  and  1.25  MeV, 
respectively)  would  be  negligibly  attenuated  [39,  15].  Prettyman  et  al.  incorporated 
information  from  a  traditional,  transmission  imaging  setup,  combined  synergistically 
with  a  Compton  scatter  tomography  setup  [56].  Such  a  technique  would  certainly 
be  viable  in  medical  imaging  applications  where  attenuation  information  is  already 
be  known.  However,  a  standalone  system  could  be  required,  especially  in  industrial 
applications  where  access  to  only  one  side  of  an  object  is  given.  With  such  constraints 
iterative  techniques  must  be  used,  such  as  those  which  were  first  rigorously  developed 
by  Hussein  et  al.  and  Arendtsz  [36,  7].  Regularization  methods  are  used  to  solve 
the  problem,  such  as  the  ML-EM  and  Penalized  Weighted  Least  Squares  (PWLS) 
methods.  Evans  utilized  the  PWLS  method  in  his  work  and  developed  a  Fortran 
code  for  its  execution  [25]. 

A  group  from  France  is  showing  some  interest  in  Compton  scatter  tomography  for 
medical  applications  [51,  52,  23,  22],  They  expand  upon  a  method  called  the  Com¬ 
pounded  Conical  Radon  Transform,  a  form  of  filtered  backprojection.  This  methodol¬ 
ogy  may  be  useful  in  further  AFIT  work  on  MCST.  Unfortunately,  they  assume  that 
attenuation  in  the  medium  is  negligible.  Earlier,  Arendtsz  argued  that  attenuation  in 
the  scattering  medium  is  the  greatest  stumbling  block  to  MCST  [7] .  Nevertheless,  the 
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Compounded  Conical  Radon  Transform  method  is  worth  note.  They  examine  some 
iterative  methods,  as  well.  They  simulate  both  collimated  and  uncollimated  MCST 
systems  in  their  work.  Their  work  is  done  with  the  Single  Photon  Emission  Computed 
Tomography  (SPECT)  application  in  mind,  making  theirs  a  subtly  different  problem 
than  the  MCST  problem  investigated  here. 

B.  Guerin  and  G.  El  Fakhri  of  the  Massachusetts  General  Hospital  and  the  Har¬ 
vard  Medical  School,  respectively,  argue  that  energy  data  already  present  in  PET 
scans  can  be  used  for  Compton  scatter  compensation  [28].  Specifically,  the  Cerium- 
doped  Lutetium  Yttrium  Orthosilicate  (LYSO)  crystals  used  in  PET  have  better 
energy  resolution  than  the  formerly  used  Bismuth  Germinate  (BGO)  crystals.  The 
11.5%  resolution  of  the  LYSO  crystals  would  not  be  adequate  for  good  MCST  spatial 
resolution  (See  Figure  19).  This  group  uses  the  data  to  correct  for  Compton  scatter 
corruption  in  the  PET  data.  They  do  not  use  the  data  in  a  standalone  manner,  but 
rather  to  correct  blur  in  the  standard  PET  image. 

An  ideal  Compton  scatter  imaging  device  would  consist  of  an  uncollimated  source 
that  shines  photons  upon  the  entire  object  and  an  uncollimated  detector  array  which 
views  the  entire  object  at  once.  Such  a  problem  would  be  highly  multiplexed.  The 
signals  from  photons  that  are  scattered  in  all  the  voxels  of  the  object  would  be  con¬ 
volved  together.  As  the  problem  becomes  more  multiplexed,  numerical  instabilities 
cause  the  reconstruction  method  to  fail.  The  multiplexed  problem  is  ill-conditioned. 

Commercial  devices  are  in  use  that  perform  Compton  scatter  imaging.  One  such 
device,  called  COMSCAN,  uses  an  X-ray  generator  and  measures  the  backscattered 
photons  [66].  This  device  is  able  to  overcome  the  mathematical  issues  of  the  Compton 
scatter  problem  by  highly  collimating  both  the  source  and  the  detectors,  as  seen  in 
Figure  8.  The  device  physically  moves,  examining  one  voxel  in  the  object  at  a  time. 
In  2010,  this  particular  device  was  successfully  used  to  analyze  historic  artifacts, 
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Figure  8.  A  diagram  of  the  COMSCAN  system,  a  commercially  available  device  which 
uses  Compton  scatter  X-rays  [66]. 

including  a  fresco,  and  Egyptian  mummy,  and  a  medieval  clasp  [30].  It  was  found 
that  Compton  scatter  imaging  was  useful  for  imaging  objects  of  large  size  and/or  high 
density.  In  either  of  those  cases,  too  few  X-rays  would  transmit  through  the  object 
in  order  to  use  conventional  X-ray  imaging  techniques. 

Compton  scatter  imaging  has  been  used  for  security  applications,  to  scan  personnel 
and  equipment  at  airports,  ports,  etc.  The  full  body,  X-ray  backscatter  scanners  used 
at  airports  (infamous  for  their  ability  to  see  through  clothes)  use  low  energy  X-rays 
that  can  barely  penetrate  the  body.  This  system  is  also  highly  collimated,  like  the 
COMSCAN  system  [8]. 

The  highly  collimated  method  is  not  ideal  for  the  medical  application,  because  the 
dose  given  to  the  patient  could  be  unsatisfactorily  high.  Unlike  with  the  airport  full- 
body  scanner,  medical  X-rays  would  have  to  be  of  a  high  enough  energy  to  penetrate 
the  body.  A  slow  scan  through  the  body  with  high  energy  X-rays  would  result  in 
too  high  of  a  dose.  Thus,  multiplexing  and  lesser  collimation  would  be  useful.  In 
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the  industrial  application  the  current  systems  are  usable,  but  if  multiplexing  were 
incorporated  then  scan-time  would  decrease  dramatically. 

2.6  Other  Applications 

2.6.1  Unmanned  robots  or  UAVs 

Unmanned  Aerial  Vehicles  (UAV)  and  unmanned  robots  are  becoming  increasingly 
popular  in  the  DoD.  The  idea  of  placing  a  position-sensitive  7-ray  detector  on  one  of 
these  vehicles  is  a  promising  one.  Already,  Sodium  Iodide  (Nal)  scintillator  crystals 
have  been  placed  aboard  small  aircraft  in  order  to  detect  uranium,  thorium,  and 
potassium  over  a  wide  area  of  land  [43] .  Position  sensitivity  is  achieved  by  flying  over 
a  large  area.  If  the  detector  aboard  the  aircraft  was  position-sensitive,  then  required 
flight  times  and  distances  would  be  reduced.  A  similar  application  would  be  to  place 
a  position-sensitive  detector  on  a  ground-based  robot. 

These  types  of  vehicles  are  not  likely  to  help  in  preventing  a  nuclear  attack.  Nei¬ 
ther  uranium  nor  plutonium  is  very  active,  meaning  that  there  would  be  few  7-rays 
to  detect,  especially  from  long  distances.  Furthermore,  a  nuclear  weapon  would  prob¬ 
ably  be  encased  in  lead  or  tungsten  shielding  to  prevent  its  detection  by  7-rays.  The 
best  application  for  unmanned  vehicles  is  the  detection  of  radioactive  material  after 
a  nuclear  reactor  accident  or  a  nuclear  weapon  detonation  [70].  Whether  by  a  nuclear 
detonation  or  a  nuclear  reactor  accident,  highly  radioactive  substances  could  be  dis¬ 
pelled  over  a  large  area  of  land.  Unmanned  vehicles  can  easily  traverse  radioactive 
areas,  sending  back  readings.  These  readings  would  give  valuable  data  to  decision- 
makers  regarding  evacuations  and  movement  of  responders  and  military  personnel. 
Also,  collection  of  radioactive  effluents  is  an  important  component  of  several  treaty¬ 
monitoring  organizations.  A  position-sensitive  system  would  help  to  guide  members  of 
these  organizations  toward  the  locations  where  radioactive  effluent  content  is  highest. 
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Lastly,  a  position-sensitive  detector  would  be  useful  inside  of  a  running  reactor 
facility  in  order  to  detect  leaks. 

2.6.2  Space  Surveillance 

Since  the  1960s,  the  United  States  has  launched  numerous  satellites  with  the 
ability  to  detect  7-rays  emitted  from  terrestrial  sources  [73].  These  satellites  are  par¬ 
ticularly  well-suited  to  detect  the  7-rays  from  nuclear  detonations.  The  satellites  are 
used  to  ensure  that  other  nations  are  adhering  to  their  treaty  obligations,  especially 
with  regard  to  nuclear  weapon  tests.  Better  position-sensitive  7-ray  detectors  could 
locate  the  source  of  the  7-rays  with  better  resolution. 

2.6.3  Space-based  Gamma-Ray  Astronomy 

Worldwide,  over  eighty  satellites  have  been  launched  with  7-ray  detection  capa¬ 
bilities.  A  good  list  of  all  the  satellites  that  have  been  used  for  7-ray  astronomy  can 
be  found  at  the  NASA  website  [46].  7-ray  detectors  are  used  to  investigate  solar 
flares  and  supernova  events  [31].  Perhaps  the  most  significant  use,  currently,  is  the 
detection  of  Gamma- Ray  Bursts  (GRB).  GRBs  are  the  most  intense  electromagnetic 
events  in  the  universe.  There  is  still  not  a  full  consensus  as  to  the  origin  of  GRBs, 
though  it  is  generally  thought  that  GRBs  are  caused  by  supernovae  or  neutron  star 
creation  events  in  distant  galaxies  [72],  GRBs  last  from  between  less  than  a  second 
to  several  minutes.  Typical  durations  are  from  20-40  seconds. 

Current  satellites  systems  used  in  7-ray  imaging  can  be  large,  heavy,  and  ex¬ 
pensive.  The  Compton  Gamma- Ray  Observatory  (CGRO),  for  instance,  weighed  a 
whopping  17  tons.  Perhaps  a  better  approach  to  7-ray  astronomy  would  be  to  deploy 
a  greater  number  of  small,  inexpensive  satellites.  Cubesat  and  Nanosat  technology  is 
becoming  increasingly  popular  [18] .  Position-sensitive  7-ray  detectors  could  be  placed 
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on  a  large  number  of  small  satellites.  These  satellites  could  be  used  to  observe  a  large 
area  of  the  sky.  The  satellites  could  also  be  used  in  a  large  array  configuration,  all 
of  them  pointed  with  the  same  field-of-view,  in  order  to  increase  image  resolution  of 
stellar  y-ray  sources. 

Coded-apertures  have  successfully  been  used  on  satellites  for  y-ray  astronomy  [60, 
54,  49].  One  possible  way  to  increase  the  views  given  to  the  satellite,  thus  decreasing 
inherent  noise,  is  to  spin  the  satellite  [32,  1],  Spinning  the  satellite  introduces  several 
complexities  in  satellite  control,  communications,  and  solar  power  acquisition. 

2.6.4  Medical  Nuclear  Imaging 

Many  researchers  have  investigated  coded-apertures  for  medical  applications.  Nu¬ 
merous  nuclear  medical  spheres  have  been  considered,  including  Single  Photon  Emis¬ 
sion  Tomography  (SPECT),  Positron  Emission  Tomography  (PET)  and  Computed 
Tomography  (CT)  [5,  33,  47,  65,  42,  14].  Often  the  thyroid  phantom  is  used  as  a 
gauge  for  the  effectiveness  of  a  given  system  [3].  Small  animals  and  capillary  tubes 
have  also  been  imaged  [4,  50].  Methods  for  decoding  include  correlation  methods, 
ML-EM,  and  PWLS.  Coded- apertures  are  most  useful  in  cases  where  a  y-emitting 
radionuclide  (typcially  Tc-99m)  is  injected  into  the  body  or  ingested  by  the  patient. 
Prominent  issues  in  medical  imaging  are  the  limited  doses  that  can  be  administered  to 
patients,  limitations  on  how  close  a  detector  can  be  placed  to  source,  and  movement 
of  the  patient  over  time.  Coded-apertures  are  less  successful  in  the  CT  case,  where 
the  source  and  the  detector  are  placed  on  opposite  sides  of  the  body,  and  varying  spa¬ 
tial  attenuation  in  the  patient’s  body  causes  spatially  varying  strength  of  the  X-rays. 
Different  methods  have  been  employed  to  account  for  vignetting  effects. 
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2.6.5  Position  Sensitivity  from  a  Single  Detector  Crystal 

There  are  a  number  of  dynamic  collimators  that  have  been  considered  in  the  past. 
Some  examples,  in  addition  to  coded-apertures,  are  the  fresnel  zone  plate,  the  rotating 
slit  collimator,  the  stochastic  aperture,  and  the  Fourier  aperture  [9] .  These  apertures 
can  be  modulated  in  time  in  order  to  gain  different  views  of  the  source. 

A  rotating  mask  collimator  has  been  successfully  mated  with  a  Nal  detector,  giving 
position  sensitivity  to  the  single  crystal  [74],  The  time-dependent  signal  given  by  the 
detector  is  compared  to  the  known,  time-dependent  modulation  of  the  collimator. 
Since  the  problem  is  constrained  in  a  known,  repetitive  way,  information  about  the 
location  and  strength  of  the  source  can  be  deduced. 


III.  Computer  Modeling 


3.1  Coded  Aperture  Ray  Tracing  Forward  Model  for  ML-EM 

A  model  was  created  in  Matlab®  to  calculate  the  coded-aperture  forward  prob¬ 
lem.  The  term  forward  problem  is  used  to  describe  the  situation  in  which  the  location, 
size,  activity,  and  other  relevant  parameters  of  the  source  plane  are  known  and  the 
detector  response  is  measured.  The  forward  model  process  is  the  opposite  of  the 
image  reconstruction  process  (or  inverse  problem)  in  which  the  detector  response  is 
known  and  the  source  plane  must  be  characterized.  The  forward  model  helps  in  image 
reconstruction. 

This  model  utilizes  the  coordinates  shown  in  Figure  15.  First,  three  relevant 
coordinates  are  specified,  being  the  location  of  a  point  source,  the  location  of  the 
programmable  collimator  and  the  location  of  the  detector  (which  is  set  as  the  origin). 
The  mask  pattern  is  input,  from  which  the  collimator  is  constructed  virtually,  in  three 
dimensions.  Two  planes  are  defined,  being  the  front  and  back  of  the  collimator.  On 
each  plane,  a  10  x  10  grid  is  carved  out,  representing  the  fronts  and  backs  of  the 
chambers  of  the  collimator.  The  size  of  the  detector  crystal  and  the  size  of  the  source 
plane  also  can  be  set.  The  resolution  of  the  source  and  the  detector  planes  are  set 
by  the  user.  In  other  words,  a  factor  is  placed  in  the  code  which  changes  how  finely 
or  coarsely  the  two  planes  are  discretized,  which  can  be  adjusted  based  on  memory 
and/or  time  constraints. 

The  forward  model’s  primary  function  is  the  calculation  of  the  matrix  A  for  use 
with  the  ML-EM  algorithm,  as  described  in  section  2.4.  Within  the  script  the  user 
defines  the  activity  of  a  point  source.  All  of  the  pixels  in  the  source  plane  and  the  de¬ 
tector  plane  are  labeled.  The  program  then  systematically  moves  through  the  source 
plane,  placing  the  source  at  each  spot  in  the  source  plane.  At  each  source  position,  the 
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Figure  9.  Two  images  of  the  detector  plane  from  the  model,  given  two  different  point 
source  positions.  Notice  how  the  features  are  thinner  in  the  right  image.  This  is  due 
to  collimator  depth,  or  vignetting,  effects,  and  is  accounted  for  in  the  model. 


program  loops  through  the  detector  pixels.  For  each  source  and  detector  pixel  combi¬ 
nation,  the  program  computes  the  vector  between  the  centers  of  the  two.  By  simple 
geometry,  the  points  of  intersection  of  this  vector  and  the  front  and  rear  collimator 
planes  are  found.  From  these  coordinates  in  the  collimator  planes,  the  chambers  tra¬ 
versed  by  the  vector  are  found.  If  the  vector  passes  through  open  chambers  of  the 
collimator  only,  then  the  number  of  counts  recorded  into  A  is  given  as: 

Aij  =  (13) 


Pa  =  exp(--pt)  (14) 

P 

where  /  is  the  number  of  photons  emitted  by  the  point  source,  d  is  the  distance 
between  detector  pixel  i  and  source  location  j .  When  the  program  checks  each  vector, 
it  records  the  number  of  grid  walls  that  the  vector  passes  through,  labeled  here  as 
N.  If  a  photon  were  to  pass  through  multiple  open  collimator  elements,  it  would 
experience  some  attenuation  in  the  steel  grid  that  separates  the  collimator  elements. 
The  fractional  attenuation  due  to  one  of  these  interior  grid  walls  is  given  as  Pa,  where 
^  is  the  mass  attenuation  coefficient  (taken  to  be  that  of  steel,  given  by  the  NIST 


30 


website),  p  is  the  density  of  the  steel  and  t  is  the  thickness  of  the  steel  [37,  34],  For 
simplicity,  t  was  set  to  a  constant,  approximate  value. 

If  the  program  finds  that  the  vector  between  a  given  source  location  and  detector 
pixel  passes  through  a  closed  collimator  element,  it  sets  the  corresponding  element  of 
A  to  zero.  All  photons  that  enter  the  closed  chambers  are  assumed  to  be  attenuated 
(See  Section  4.1  for  justification  of  this  assumption).  It  was  found  experimentally  that 
a  significant  number  of  photons  are  Compton  scattered  inside  of  the  collimator.  These 
events  are  ignored  by  the  model  because  an  energy  discriminating  detector  was  used. 
The  detector  can  be  set  to  record  only  those  detected  photons  that  were  received 
with  the  known  energy  of  the  monoenergetic  radionuclide  source  (plus  or  minus  2-3 
keV).  Since  the  detector  filters  out  the  Compton  scattered  photons,  the  model  need 
not  include  them.  Additionally,  every  tenth  row  and  column  of  the  detector  matrix  is 
set  to  zero,  as  well  as  the  corner  pixels,  to  mimic  the  real  detector  (see  Section  4.2). 


Figure  10.  An  illustration  of  the  vignetting  effect.  Viewed  from  the  source,  through 
the  collimator.  A  chamber  directly  in  line  with  the  source  and  a  given  detector  pixel 
allows  full  throughput.  Off-axis  chambers  block  part  of  the  possible  paths,  as  a  result 
of  collimator  depth. 
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One  bit  of  functionality  added  later  to  the  model  was  the  ability  to  characterize, 
more  accurately,  the  vignetting  that  occurs  from  sources  in  the  near-held  region  due 
to  the  collimator’s  physical  depth.  Some  work  has  been  done  to  account  for  the  effect 
of  collimator  depth  on  correlation  methods,  particularly  by  Mu  and  Liu  in  2006  [50] . 
Using  a  different  method,  this  model  accounts  for  vignetting  by  first  over-resolving  the 
detector  plane.  The  script  allows  for  the  user  to  define  how  fine  the  resolution  should 
be.  Each  detector  plane  pixel  is  thus  split  into  smaller  squares.  The  probability  of  a 
photon  reaching  each  smaller  square  is  found  in  the  same  manner  as  before.  All  of 
the  probabilities  of  the  little  squares  are  then  averaged  to  find  the  overall  probability 
for  a  given  detector  plane  pixel. 


Table  1.  Run-times  and  memory  requirements  for  the  calculation  of  the  A  matrix. 
Run-times  are  based  on  an  above-average  desktop  computer  (for  2011). 


M 

N 

Time 

RAM  Req.  (MB) 

100 

25,600 

10m 

10 

2,500 

25,600 

4h7m 

250 

3,600 

25,600 

6h 

369 

10,000 

25,600 

16h 

1,024 

The  computation  of  the  A  matrix  is  nontrivial.  As  M  increases,  the  RAM  and 
time  required  to  compute  A  also  increases,  as  seen  in  Table  1.  There  is  a  limit  to 
how  finely  the  source  and  detector  planes  can  be  discretized  based  on  the  computer 
setup. 

In  a  future  version  of  the  model,  differences  in  strip  efficiencies  should  be  taken  into 
account,  as  well  as  the  change  in  efficiency  across  a  strip.  As  with  the  MOST  forward 
model,  inclusion  of  detector  depth  information  would  increase  model  accuracy. 
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3.2  MCNP  Forward  Model 


A  forward  model  was  also  created  using  a  Monte  Carlo  program  called  Monte 
Carlo  N-Particle  Transport  Code  (MCNP).  MCNP  was  originally  developed  by  Los 
Alamos  National  Laboratory  to  model  neutron  transport  in  nuclear  reactors,  but  it 
can  be  useful  for  modeling  7-rays  and  their  interactions  with  matter  as  well  [69]. 
MCNP  was  chosen  over  other  Monte  Carlo  programs  simply  because  of  the  author’s 
familiarity  with  the  code  and  its  immediate  availability.  A  primer  by  Shultis  and  Faw 
of  KSU  was  frequently  referenced  during  the  writing  of  this  model  [63] . 

The  term  Monte  Carlo  was  first  coined  and  used  extensively  by  John  von  Neu¬ 
mann,  Stanislaw  Ulam,  and  Nicholas  Metropolis  for  the  Manhattan  Project  [24],  In 
Monte  Carlo  methods,  the  geometry  of  the  system  is  carefully  specified  so  that  the 
material  of  all  of  the  objects,  or  cells,  in  the  system  is  specified.  Starting  at  the  source 
position,  dice  are  rolled,  figuratively,  and  a  photon  is  emitted  in  a  random  direction. 
A  large  number  of  distributions  are  stored  in  the  software  for  many  different  types 
of  materials.  The  distributions  are  based  on  cross  sections  for  interactions.  For  in¬ 
stance,  the  software  has  a  stored  distribution  describing  how  often  a  Compton  scatter 
event  will  occur  in  air  for  a  photon  of  a  given  energy.  The  code  will  randomly  pick 
different  outcomes  for  each  photon,  weighted  by  the  appropriate  distributions.  An 
outer  sphere  is  usually  defined,  the  interior  of  which  is  called  the  universe.  Photons 
that  exit  the  sphere  are  “killed.”  The  code  probabilistically  determines  the  path  of 
a  single  photon  until  that  photon  is  either  killed  or  deposits  all  of  its  energy  in  one 
of  the  cells  within  the  universe.  Once  either  of  those  criteria  are  met,  the  program 
records  the  path  of  the  photon  and  starts  a  new  photon.  Millions  of  photons  can  be 
tracked  in  this  fashion. 

A  Matlab®  script  was  developed  to  write  the  input  hie  for  MCNP.  The  in¬ 
put  hies  were  thousands  of  lines  long  and  therefore  could  not  reasonably  be  written 
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Figure  11.  Left:  The  programmable  collimator,  as  modeled  in  MCNP.  Right:  A  side 
view  of  the  programmable  collimator  and  the  detector  crystal,  as  modeled  in  MCNP. 

manually.  To  use  the  Matlab®  script,  the  user  simply  inputs  the  location  of  the 
source  and  the  mask  relative  to  the  detector  crystal  in  three  dimensions,  as  well  as 
the  mask  configuration.  Based  on  the  mask  configuration,  the  script  fills  in  the  cells 
corresponding  to  the  collimator  chambers  either  with  air  and  a  PEEK  plug  at  the 
rear,  or  with  a  PEEK  plug  at  the  front  and  AIM-70  in  the  rest  of  the  chamber, 
depending  on  whether  the  chamber  is  “on”  or  “off.”  The  script  also  defines  160  x 
160  detector  cells,  corresponding  to  the  detector  subpixels.  A  122  keV  y-ray  point 
source  was  used.  Initially  the  detector  crystal  was  made  from  germanium.  The  idea 
was  to  record  tallys.  Each  time  a  photon  deposits  its  full  energy  in  a  detector  cell, 
a  tally  corresponding  to  that  cell  increments  up  by  one.  However,  MCNP  only  will 
track  1000  tallys,  significantly  less  than  the  25,600  cells  of  the  detector  crystal.  After 
this  realization,  it  was  decided  that  the  ptrac  hies  would  be  used.  The  ptrac  hies  list 
the  track  of  all  of  the  particles  of  interest.  The  detector  crystal  material  was  set  to 
vacuum,  and  any  photon  that  passed  through  the  detector  in  the  appropriate  energy 
range  is  recorded  in  the  ptrac  hie.  Just  like  with  the  real  detector,  MCNP  has  an 
energy  window,  which  was  set  to  122  keV  +/-  5  keV.  At  that  time  the  plan  was  that 
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the  actual  detector  would  use  a  10  keV  energy  window,  though  this  was  later  reduced 
to  6  keV  to  further  filter  out  Compton  scattered  photons. 

To  further  increase  computational  efficiency,  the  photons  emitted  from  the  source 
were  constrained  to  a  cone  which  subtended  an  area  slightly  bigger  than  the  pro¬ 
grammable  collimator.  In  other  words,  the  source  was  only  allowed  to  emit  photons 
in  the  solid  angle  of  the  cone.  By  not  allowing  photons  to  be  emitted  in  the  oppo¬ 
site  direction  of  the  programmable  collimator,  for  instance,  computation  time  was 
reduced  dramatically. 


Figure  12.  Predictions  of  the  detector  response  from  a  point  source,  given  by  MCNP, 
using  100,000  particles  (Left)  and  1,000,000  particles  (Right). 


A  ptrac  parser  was  written  to  interpret  the  ptrac  data.  For  a  given  particle  track, 
the  parser  records  which  detector  cells  that  the  path  traversed.  Also,  for  a  given  path, 
the  parser  records  the  number  of  detector  cells  traversed,  called  EC  (for  event  count). 
The  parser  script  then  adds  to  each  of  the  affected  detector  cells,  assuming  full 
energy  deposition  to  be  equally  possible  for  each  of  the  traversed  cells.  This  is  done 
for  all  particle  tracks. 

At  one  time,  there  was  hope  that  MCNP  could  be  used  to  provide  the  A  matrix 
for  ML-EM.  However,  without  a  sufficient  number  of  particles  in  MCNP,  there  is  too 
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Table  2.  Run-times  and  ptrac  file-sizes  for  MCNP  simulations  based  on  the  number  of 
particles  emitted  from  the  source.  Run-times  are  based  on  an  above-average  desktop 
computer  (for  2011). 


Particles 

File-Size  (MB) 

Run-Time 

Avg.  Counts  per  Det.  Subpixel 

10,000 

0.8 

13.6s 

0.00377 

50,000 

2.7 

43s 

0.018 

100,000 

5 

lml8s 

0.0329 

500,000 

24 

5m  16s 

0.1736 

1,000,000 

47 

9m57s 

0.349 

5,000,000 

235 

47ml6s 

1.74 

10,000,000 

471 

lh33m 

- 

100,000,000 

4,698 

15h22m 

- 

much  variation  from  pixel  to  pixel  in  A.  The  data  must  be  smooth  if  they  are  to  work 
for  ML-EM.  To  achieve  the  proper  smoothness,  at  least  10,000,000  particles  must  be 
used,  although  100,000,000  would  be  better.  An  average  A  matrix  has  about  1,000 
source  positions  with  up  to  10,000  in  some  cases.  Conservatively,  using  the  data  from 
Table  2,  it  would  take  65  days  to  compute  A.  This  simply  is  not  reasonable,  especially 
given  the  fact  that  a  different  A  is  needed  for  every  source-collimator-detector  ge¬ 
ometry  and  mask  configuration.  Nevertheless,  the  MCNP  simulation  provided  good 
insight  into  the  number  of  expected  Compton  scatter  events,  the  expected  attenua¬ 
tion  through  the  “off”  chambers,  and  things  along  similar  lines.  Perhaps  by  using 
variance  reduction,  or  other  more  sophisticated  techniques,  MCNP  could  be  used  to 
calculate  the  forward  model  for  ML-EM. 

3.3  The  Optimal  Mask  Sequence 

A  code  was  developed  to  find  the  optimal  mask  sequence  for  the  coded-aperture 
application  of  the  programmable  collimator.  The  code  iteratively  determines  a  se¬ 
quence  of  matrices  that  are  most  independent  from  each  other.  The  cross-correlation 
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Create  the  initial  sequence  of  random 
mask  matrices 


I 


v  random  mask 


Place  the  mask  at  each  position  in  the 
sequence 


Find  the  total  correlation  (the  sum  of  the 
correlation  between  each  mask  and  all 
other  masks) 


Determine  which  mask  position  caused 
the  lowest  total  correlation. 


Place  mask  into  sequence  at  position  found 
above,  if  the  total  correlation  of  the  sequence 
is  less  with  the  new  mask  than  without. 
Otherwise,  sequence  stays  unchanged. 


Create  a  sequence  of  random  mask 
matrices 


Find  the  total  correlation  (the  sum  of  the 
correlation  between  each  mask  and  all 
other  masks) 


If  the  total  correlation  is  lower  than  that  of 
the  previous  iterations  sequence,  store  the 
new  sequence.  Otherwise  continue  to  use 
the  old  sequence 


Figure  13.  Methods  to  find  the  sequence  of  masks  that  are  least  correlated  with  each 
other. 


is  used  as  a  measure  of  independence  between  mask  matrices.  More  specifically,  the 
magnitude  of  the  central  spike  of  the  cross-correlation  is  used  as  the  metric.  If  the 
central  spike  is  zero,  then  the  two  mask  matrices  are  completely  independent.  The 
two  iterative  methods  are  shown  in  Figure  13. 

It  was  found  that  the  right-hand  method  shown  in  Figure  13  converged  more 
quickly  than  did  the  left-hand  method.  Also,  by  observing  the  operation  of  these 
codes,  it  was  discovered  that  the  least  correlated  sequence  of  matrices  is  that  in 
which  each  mask  clement  is  open  or  closed  for  equal  durations.  For  instance,  in  a 
sequence  of  six  mask  configurations,  each  element  should  be  open  in  three  of  the 
masks  and  closed  in  three  of  the  masks.  After  this  discovery  was  made,  a  code  was 
written  to  construct  a  sequence  of  random  masks  that  adhere  to  this  criterion. 

An  additional  algorithm  was  written  which  orders  the  masks  in  the  sequence  for 
which  the  programmable  collimator  user  will  have  the  least  number  of  element  changes 
to  make  during  an  experiment.  This  algorithm  can  lower  the  number  of  changes  by  40- 
50  for  a  six-mask  sequence  of  10x10  matrices.  Later,  code  may  be  developed  in  which 
the  number  of  changes  is  determined  by  the  user,  and  the  code  will  find  the  sequence 
of  least-correlated  masks  with  the  specified  number  of  changes.  Additionally,  this 
code  could  be  modified  to  choose  the  best  masks  on  the  fly.  Incoming  information 
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would  be  incorporated  to  help  choose  the  next  mask  to  use.  These  optimal  mask 
sequences  would  vary  based  on  the  application. 


Figure  14.  An  example,  optimized  sequence  of  random  mask  matrices.  Notice  that 
each  mask  element  in  this  six-mask  sequence  is  open  in  three  masks  and  closed  in  three 
masks. 


3.4  MCST  Forward  Model 

A  three  dimensional  Cartesian  coordinate  system  was  devised.  The  millimeter 
was  chosen  as  the  unit  in  this  space.  The  detector  crystal  is  assumed  to  be  a  square, 
and  the  origin  was  placed  at  the  furthest,  lowest  point  of  the  detector  to  keep  all 
coordinates  positive.  Only  three  points  in  space  need  to  be  defined  for  a  given  photon 
event:  the  point  source,  a  point  in  the  material,  and  a  point  in  the  detector.  Vectors 
following  the  rays  of  light  are  found  between  the  source  and  the  material  and  between 
the  material  and  the  detector,  defined  as 

srcmatvect  (x mat  src,  ymat  Vsrci  % mat  zSrc)  (15) 

raatpixvect  (%pix  mati  Vpix  2/mat)  % pix  % mat )  (16) 
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The  Compton  scattering  angle,  9,  is  found  through  the  dot  product,  such  that 


Of  course,  the  detector  pixel  is  not  a  point,  but  rather  is  three  dimensional.  Each 
pixel  is  approximated  as  a  square  that  extends  in  y  and  z.  The  angles  from  the 
point  material  to  the  top  and  bottom  of  the  detector  pixel,  called  /3max  and  / 3min ,  are 
found  using  the  same  methodology  as  used  above.  In  the  same  way,  the  angles  to 
the  left  and  right  sides  of  the  detector  pixel,  amax  and  am.in:  are  found.  These  angle 
ranges  are  used  to  find  the  solid  angle  subtended  by  the  detector  pixel  as  part  of  the 
Klein-Nishina  integration. 

A  point  material  was  placed  in  the  space  based  on  the  coordinates  in  Figure  15. 
A  detector  was  constructed  virtually  as  an  80x80  mm2  square,  with  256,  5x5  nun2 
pixels.  The  scattered  energy,  E\  was  calculated  using  the  Compton  equation  for  each 
detector  pixel  (Equation  12).  Given  a  Z  of  29  for  copper,  the  expected  cross  section, 
in  barns,  for  each  detector  pixel  was  calculated. 


Looking  from  Door  Side  of  Dome 


Figure  15.  The  coordinate  system  and  relevant  dimensions  used  for  modeling  with  the 


PHDs  Co.  DSSD 
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Eventually,  this  model  can  be  used  in  conjunction  with  an  ML-EM  or  PWLS 
method.  It  may  be  useful  to  include  detector  depth  information  in  this  model.  The 
PHDs  strip  detector  system  is  able  to  make  relatively  accurate  estimates  of  the  depth 
at  which  an  event  occurred  [75].  The  inclusion  of  depth  information  would  result  in 
a  more  accurate  model,  geometrically.  And,  of  course,  the  more  accurate  the  model 
is,  the  better  the  iterative  methods  like  ML-EM  and  PWLS  will  perform. 


Model  Prediction  for  energy  of  y  scattered  from  point  Cu  material.  Colorbar  in  keV. 
1 1 1018_MCST_Cd109_CuStrip  geometry  used. 


2  4  6  8  10  12  14  16 


Detector  Column 


Model  Prediction  for  cross  section  of  y  scattered  from  point  Cu  material  into  detector  pixel. 

Colorbar  in  bams.  111018_MCST_Cd109_CuStrip  geometry  used.  x  k 
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Figure  16.  Left:  The  predicted  photon  energy  into  each  pixel  from  the  copper  “point” 
phantom.  Right:  The  cross  section  for  scatter,  from  the  copper  “point”  phantom  into 
each  detector  pixel. 


3.5  MCST  Resolution  and  Setup  Analysis 

At  first,  when  thought  was  given  to  an  MCST  setup,  it  was  unclear  what  the 
optimal  setup  would  be.  On  the  one  hand,  the  Klein-Nishina  distribution  is  angle 
dependent.  The  scattering  cross-section,  in  the  backscatter  region,  is  greatest  at  180 
degrees.  So  to  get  the  greatest  number  of  scatters  over  a  given  time,  the  source  and 
the  detector  would  have  to  be  close  to  each  other  with  the  sample  between  them 
in  y  and  out  a  ways  in  x.  On  the  other  hand,  as  the  sample  moves  closer  to  the 
detector  the  angles  of  the  source-detector  pixel  vectors  are  differentiated  more.  With 
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Figure  17.  Left:  The  energy  of  an  88  keV  photon  after  a  Compton  scatter,  from  the 
Compton  formula.  Right:  The  differential  cross  section  for  scatter  in  various  materials, 
from  the  Klein-Nishina  formula  [37]. 


greater  differences  in  the  angles  into  each  detector  pixel,  there  are  greater  differences 
in  the  energies  found  in  each  detector  pixel.  A  goal  of  the  detector  positioning  is 
that  the  detector  subtends  the  greatest  solid  angle  of  scattered  photon  paths  from 
the  material. 


Figure  18.  At  the  closer  material-detector  distance,  X\,  the  detector  subtends  a  much 
larger  solid  angle,  cf>  as  compared  to  the  solid  angle  6  subtended  by  the  detector  when 
the  material-detector  distance  is  increased  only  modestly  ( x 2). 
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It  was  discovered  empirically  that  the  angle  discrimination  advantage  gained  by 
moving  the  material  closer  to  the  detector  outweighs  the  decreased  cross  section 
for  scatter.  Indeed,  to  allow  for  the  closest  material-detector  distance,  the  source- 
material-detector  angle  becomes  about  90  degrees,  being  the  angle  with  the  lowest 
scattering  cross  section.  In  a  real  application,  the  sample-detector  distance  may  be 
fixed  to  a  longer  than  desirable  length.  To  get  the  best  results  this  distance  was 
minimized  in  experiments.  The  effect  of  this  distance  on  potential  resolution  is  given 
in  Figure  19. 


The  spatial  resolution  possible  given  an  average  detector  response  energy  difference, 


For  different  detector-material  distances.  Based  on  the  Compton  Fomula 


The  difference  in  energy,  at  an  average  pixel,  between  the  two  detected  Gaussian  peaks  [keV] 


Figure  19.  Using  the  model,  two  infinitesimally  small  pieces  of  copper  were  separated 
by  various  distances  (y-axis).  The  difference  in  Compton  scattered  photon  energy  from 
each  piece  was  found  (x-axis).  This  graph  gives  the  MCST  spatial  resolution  possible 
in  the  sample  given  the  minimum  difference  in  energy  between  detected  photons  that 
is  distinguishable.  Three  sample-detector  distances  are  shown. 


To  make  things  more  clear,  a  two-dimensional  phantom  was  conceived,  as  shown 
in  Figure  20.  A  model  was  written  to  predict  the  Gaussian  response  that  each  piece 
of  material  would  cause  in  the  detector.  These  responses  must  be  distinguishable  if 
the  MCST  method  is  to  be  effective.  As  seen  in  Figure  20,  the  distance  between 
the  material  and  the  detector  significantly  affects  this  distinguishability.  To  reiterate, 
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the  minimization  of  the  detector-sample  distance  should  be  the  primary  focus  when 
considering  the  dimensions  of  the  MOST  setup.  In  importance,  the  minimization  of 
this  distance  far  outweighs  the  maximization  of  the  cross  section  for  scatter. 


Figure  20.  Top:  The  AlCuCu  phantom  used  for  modeling.  Left:  The  predicted  detector 
response,  if  the  phantom  were  located  at  105.4  mm.  Right:  The  predicted  detector 
response  if  the  phantom  were  located  at  205  mm.  Notice  how  indistinguishable  the 
Gaussians  from  the  three  materials  become  when  the  detector /material  distance  is 
doubled. 
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IV.  Design  and  Experimental  Setup 


4.1  Progammable  Collimator  Design 

In  the  first  design  iteration,  it  was  thought  that  plastic  rods  would  be  submerged 
in  a  liquid  metal,  all  within  a  sealed  container.  The  rods  would  have  small  magnets  on 
each  end  and  would  be  actuated  by  external  permanent  magnets  or  electromagnets. 
However,  the  attenuation  in  the  plastic  rods  was  calculated  as  being  too  great  and 
significant  difficulties  were  foreseen  in  actuating  the  rods.  Eventually  the  current 
design  was  chosen. 


Figure  21.  Molten  AIM-70 


The  programmable  collimator,  in  its  present  form,  was  built  by  the  AFIT  model 
shop  in  the  Summer  and  Fall  of  2011.  A  steel  block  was  fabricated,  using  Electrical 
Discharge  Machining  (EDM),  to  include  100  elements,  or  chambers,  in  a  10x10  array. 
The  square,  cross-section  of  each  chamber  is  5x5  mm2,  like  the  PHDs  DSSD  detector 
pixels.  Each  chamber  can  be  turned  “on”  or  “off”  toy-rays.  The  design  was  optimized 
for  the  511  keV  y-rays  from  electron/positron  annihilation.  Originally  the  device 
was  intended  for  AFIT’s  ongoing  positron  work.  Lower  and  higher  energy  photons 
certainly  can  be  imaged  through  this  collimator. 
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Figure  22.  Left:  The  front  face  of  the  programmable  collimator.  Right:  The  tool  used 
to  move  the  PEEK  rods  back  and  forth  in  the  collimator  chambers,  locked  to  a  single 
PEEK  plug. 


The  primary  component  of  the  collimator  apparatus  is  the  grid  block,  as  shown 
in  the  left  of  Figure  23.  The  block’s  outer  dimensions  are  5x5x5  cm3.  The  block 
is  suspended  in  the  center  of  the  outer  housing.  While  in  operation  the  housing  is 
filled  with  a  liquid  metal  called  AIM-70,  also  known  as  CerroBEND,  Wood’s  Metal  or 
MCP  158.  AIM-70  is  a  eutectic  alloy  composed  of  50%  bismuth,  26.7%  lead,  13.3% 
tin  and  10%  cadmium,  by  weight,  with  a  melting  point  of  only  70°  C  (158°  F)  [48]. 
Due  to  its  high  Effective,  AIM-70  is  a  good  attenuator  of  y-rays.  Each  chamber  also 
contains  a  plug  made  of  Poly  Ether  Ether  Ketone  (PEEK),  which  prevents  the  liquid 
AIM-70  from  spilling  out  of  the  front  of  the  grid  block.  PEEK  is  a  highly  machinable 
plastic  that  is  able  to  withstand  high  temperatures.  The  particular  flavor  of  PEEK 
chosen  is  filled  with  carbon  to  30%.  By  filling  the  PEEK  with  carbon,  the  com¬ 
pressive  strength  increases,  the  expansion  rate  decreases  and  the  themal  conductivity 
increases,  as  compared  to  unfilled  PEEK  [55].  Having  a  low  %e ff active,  PEEK  is  a 
poor  7  attenuator. 

Centered  on  the  rear  plate  of  the  apparatus  is  a  square  window,  also  made  of 
PEEK.  This  window  is  5  mm  thick.  When  the  apparatus  is  assembled,  the  rear  face 
of  the  chamber  block  is  only  3  mm  from  the  rear  PEEK  window.  This  spacing  allows 
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Figure  23.  Left:  The  rear  face  of  the  chamber  block,  suspended  in  the  center  of  the 
housing  plate,  as  seen  from  the  interior.  Right:  The  PEEK  window  in  the  back  of  the 
collimator  assembly. 

the  AIM-70  to  flow  into  the  chambers  from  the  outer  housing  box.  A  “T”  pattern 
was  cut  out  of  each  PEEK  plug  (Figure  22).  A  matching  tool  was  fabricated  with  a 
similar  “T”  shape  on  the  end.  The  tool  is  inserted  into  the  chamber  with  the  proper 
orientation  to  allow  entry  into  the  cutout  of  the  PEEK  plug.  Once  inside  the  plug, 
the  tool  is  turned  ninety  degrees,  locking  it  in  place.  The  user  may  then  push  and 
pull  the  PEEK  plug  within  the  chamber. 

To  operate  the  collimator,  each  of  the  100  chambers  is  configured  in  one  of  two 
ways.  The  chamber  is  turned  “on,”  meaning  that  y-rays  are  allowed  through,  or 
“off,”  meaning  that  y-rays  are  attenuated.  To  turn  the  chamber  “off,”  the  PEEK 
plug  is  pulled  all  the  way  to  the  front  of  the  grid  block.  The  AIM-70,  which  fills 
the  outer  housing,  flows  around  the  back  of  the  grid  block  into  the  chamber.  The 
AIM-70  in  the  chamber  attenuates  any  y-rays  that  try  to  pass  through.  To  turn  a 
given  chamber  “on,”  the  PEEK  plug  is  pushed  as  far  back  as  possible  so  that  it  mates 
with  the  rear  PEEK  window.  While  the  PEEK  plug  is  pushed  backward,  the  AIM-70 
is  displaced  from  the  chamber  and  into  the  outer  housing  assembly.  Once  the  PEEK 
plug  is  mated  to  the  rear  PEEK  window,  a  y-ray  passing  through  that  particular 
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3  mm  of  space  b/t  chamber  block  and 
rear  PEEK  window  where  the  AIM-70 
can  flow  into  the  chambers 


Rear 

PEEK 

Window 


Figure  24.  A  slice  through  the  programmable  collimator.  If  the  PEEK  plug  is  pulled 
out  (top  chamber),  the  chamber  fills  with  AIM-70,  and  7-rays  are  blocked.  If  the  PEEK 
plug  is  pushed  in  (bottom  chamber),  the  AIM-70  is  pushed  out  of  the  chamber  and 
7-rays  must  only  pass  through  a  small  bit  of  PEEK. 
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chamber  will  encounter  only  air  and  17  mm  of  PEEK,  and  thus  will  be  only  lightly 
attenuated. 

More  quantitatively,  it  is  known  that 

=  e_(p)pt  (18) 

lo 

where  I  is  the  number  of  transmitted  photons,  Iq  is  the  original  number  of  photons, 
is  the  mass  attenuation  coefficient  and  pt  is  the  mass  thickness  of  the  material 
in  question  [37].  Hubbell  and  Seltzer’s  NIST  data  were  used,  paim- 70  was  taken 
as  9.38  ^3  and  (>peek  was  assumed  to  be  p Polyethylene,  which  is  1.19  [34],  As 

one  example,  at  511  keV,  83%  of  the  photons  will  pass  through  a  chamber  in  the 
“on”  configuration,  whereas  only  0.57%  of  the  photons  will  pass  through  a  chamber 
in  the  “off”  configuration.  For  100  keV  photons,  72%  of  the  7-rays  would  make  it 
through  an  “on”  chamber,  whereas  10-70%  of  the  7-rays  would  make  it  through  an 
“off”  chamber. 


Energy  (keV) 

Length  [cm] 

AIM70  (mu/rho) 
[cmA2/g] 

AIM-70  I/I  0 

PEEK(mu/rho) 

[cmA2/g] 

PEEKI/I  0 

100 

1.7 

4.~3E+00 

1.85546E  33 

0.1641 

0.71T50576 

100 

3.6 

4.73E+00 

4.85438E-70 

0.1641 

0.4950956 

511 

1.7 

0.1476 

0.095023638 

0.0941 

0.82665793 

511 

3.6 

0.1476 

0.006845548 

0.0941 

0.66822897 

1000 

1.7 

0.0686 

0.334909277 

0.0687 

0.87024535 

1000 

3.6 

0.0686 

0.098619645 

0.068- 

0.7450449 

Figure  25.  Selected  attenuation  data  for  AIM-70  and  PEEK  (approximated  as  polyethy¬ 
lene)  for  two  different  lengths  at  three  different  energies. 


To  keep  the  AIM-70  in  liquid  form,  the  entire  apparatus  must  be  kept  above 
70°  C.  This  was  accomplished  using  a  hot  plate.  In  the  future,  heating  tape  may 
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be  used  to  maintain  operating  temperatures.  Initially  there  was  concern  that  the 
PEEK  plugs  would  fit  differently  within  the  steel  chambers  upon  heating,  in  that 
the  thermal  expansion  coefficient  for  each  material  is  different.  In  other  words  PEEK 
and  steel  expand  at  different  rates  when  heated.  In  practice,  however,  this  differential 
expansion  seems  to  be  negligible.  The  plugs  move  similarly  at  all  temperatures  tested. 

An  initial  test  with  AIM-70  resulted  in  some  leakage  from  the  housing.  The  outer 
housing  was  disassembled  and  the  edges  that  make  contact  with  each  other  were 
ground  and  polished,  resulting  in  no  further  leakage.  In  another  test,  the  apparatus 
was  filled  with  water  under  the  assumption  that  water  would  leak  more  than  would 
the  thicker,  liquid  AIM-70.  The  apparatus  was  heated  to  operating  temperatures  and 
the  PEEK  plugs  were  pushed  and  pulled  within  the  chambers.  During  this  test,  no 
water  leaked  from  any  of  the  junctures.  Furthermore  this  test  again  proved  that  the 
difference  in  expansion  between  the  PEEK  and  steel,  with  temperature  increase,  is 
inconsequential. 

As  the  collimator  was  tested,  concern  developed  that  the  experimenter  might 
pull  a  given  PEEK  plug  all  the  way  out  of  the  grid  block  when  trying  to  change 
a  given  chamber  to  the  “off”  configuration.  If  this  were  to  occur,  molten  AIM-70 
would  eject  from  the  device,  spilling  out  into  the  experimental  area.  To  assuage  this 
risk,  the  collimator  configuration  was  changed  in  a  fume  hood.  The  collimator  was 
placed  within  a  large,  metal  tub.  The  AIM-70  happened  to  cool  and  solidify  by 
the  time  it  was  taken  to  the  detector  setup.  The  cooling  AIM-70  is  something  that 
warrants  consideration.  AIM-70  expands  0.006  inches  per  inch  upon  cooling  [13]. 
This  expansion  causes  stress  on  the  chamber  walls  within  the  grid  block.  A  few  of 
the  chamber  walls  have  shown  cracking,  possibly  due  to  this  effect. 

Ultimately,  a  different  alloy  is  desired.  Another  eutectic  alloy  called  AIM-47  (or 
LOW  117)  may  be  considered.  It  has  roughly  the  same  ZEffective  as  AIM-70  but  has 
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a  lower  melting  temperature  (47°  C  or  117°  F)  and  less  expansion  upon  solidification 
(+0.0002  inches  per  inch  over  the  first  6  minutes,  followed  by  shrinkage  to  -0.0002 
inches  per  inch).  This  is  discussed  further  in  Section  6.4. 

The  FOV  of  the  collimator  was  constrained  because  of  the  collimator’s  depth.  The 
height  or  width  of  the  square  plane  visible  to  the  collimator,  HknyChamb,  is  given  in  cm 
by 

2D 

-^AnyChamb  5  (19) 

where  D  is  the  distance  from  the  front  of  the  collimator  (or  the  source  side  of  the 
collimator)  to  the  source  plane,  in  cm.  The  equation  gives  the  height  of  the  maximum 
visible  plane,  meaning  that  source  positions  in  which  the  source  is  visible  through  as 
little  as  a  single  chamber  are  included.  If  the  source  is  to  be  seen  through  all  of  the 
chambers,  then  the  visible  square  is  somewhat  smaller,  where  the  width  or  height, 
^Aiichamb,  in  cm,  is  given  by 

AllChamb  AnyChamb  9  )  D  20  (20) 

Notice  that  the  source  must  be  at  least  20  cm  from  the  front  of  the  collimator  if  it  is 
to  be  viewable  through  all  of  the  chambers. 

4.2  The  PHDs  HPGe  Strip  Detector 

A  PHDs  Co.  Double-Sided  Strip  Detector  (DSSD)  was  used  in  the  experimental 
portions  of  this  study.  The  detector  consists  of  a  single,  9  cm  diameter  High-Purity 
Germanium  (HPGe)  crystal,  with  an  active  depth  of  11  mm.  The  crystal  is  split  into 
16  strips  on  the  front  and  rear  faces  of  the  detector,  where  the  front  and  rear  strips 
are  orthogonal  to  each  other.  Because  the  crystal  is  round,  only  the  middle  ten  strips 
on  each  side  extend  to  a  length  of  8  cm.  Moving  outward  from  the  center  ten  strips, 
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each  strip  is  progressively  5  mm  shorter,  so  as  to  conform  to  the  shape  of  the  crystal. 
The  strips  are  5  mm  wide,  effectively  resulting  in  220,  5x5  mm2  energy-sensitive  pixels 
(Figure  26).  The  rear  strips,  which  run  vertically,  are  termed  the  AC  side,  whereas 
the  front  strips,  which  run  horizontally,  are  termed  the  DC  side.  The  AC  strips  are 
labeled  from  0-15,  where  0  is  the  left-most  strip  and  15  is  the  right-most  strip.  The 
DC  strips  are  labeled  from  16-31,  where  16  is  the  top  strip  and  31  is  the  bottom  strip 
[75].  The  segmented  detector  contacts  were  fabricated  using  amorphous-germanium 
detector  technology  [29,  44,  57,  35].  A  Stirling-cycle  mechanical  cooler  maintains 
the  detector  operating  temperature  near  68  K  in  a  vacuum  crystat  assembly  having 
a  mass  of  about  12  kg  [78].  The  mechanical  cooling  and  low  weight  mean  that  the 
programmable  collimator  and  detector  combination  could  be  made  portable  for  future 
applications. 


I 

5  mm 
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Figure  26.  The  configuration  of  the  PHDs  detector  crystal  [75]. 

The  strips  are  connected  to  an  electronics  system  called  the  Spect32.  The  Spect32 
contains  four  Field  Programmable  Gate  Arrays  (FPGA),  each  optimized  for  imag¬ 
ing.  Subpixcl  resolution  can  be  attained  through  subpixel  interpolation.  When  the 
electrons  and  holes  move  through  a  given  strip  after  a  photon  event,  transient  charge 
is  induced  in  the  adjacent  strips.  The  amount  of  transient  charge  induced  in  one 
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adjacent  strip  relative  to  the  other  adjacent  strip  is  related  to  the  location  across  the 
event  strip  in  which  the  photon  event  occurred.  If  the  photon  event  occurred  closer 
to  one  side  of  the  event  strip,  more  transient  charge  will  be  induced  in  the  closer 
adjacent  strip  [11,  75,  67].  The  transient  charge  in  the  adjacent  strips  is  measured 
by  the  detector  system.  By  interpolation,  subpixcl  resolution  down  to  0.5  mm  is 
achieved.  Thus,  the  overall  resolution  of  the  detector  is  increased  by  a  factor  of  ten. 
Several  subpixel  interpolation  schemes  are  built  in  to  the  detector  system  software. 
Every  tenth  row  and  column  subpixcl  always  reads  zero  counts,  corresponding  to  the 
location  of  the  physical  separator  between  strips. 

The  detector  has  a  FWHM  energy  resolution  of  1  keV  at  88  keV.  Before  use  the 
detector  must  be  calibrated  in  energy  using  a  known  source.  Ideally,  the  same  source 
that  is  to  be  imaged  should  be  used  for  calibration. 

4.3  MCST  Uncollimated  Forward  Problem  Experimental  Setup 

A  view  from  the  top  of  the  basic  setup  for  MCST  is  shown  in  Figure  27.  A  101.8 
yuCi  (on  07NOV11),  Cd-109  button  source  was  used  (T-098).  The  source  was  secured 
to  a  wooden  block,  and  the  block  was  secured  to  a  precision  lab  jack.  A  set  of  three 
girders,  arranged  as  three  sides  of  a  square,  was  placed  in  front  of  the  detector.  The 
square  formed  by  the  girders  was  much  taller  and  wider  than  the  detector,  so  that 
Compton  scattered  photons  from  the  girders  would  be  attenuated  geometrically  and 
so  that  Compton  scattered  photons  from  the  girders  would  be  of  a  much  different 
energy  than  those  from  sample.  Two  pieces  of  nylon  fishing  line  then  were  hung  from 
the  top  girder.  Weights  were  tied  to  the  bottom  of  each  fishing  line,  to  keep  the 
lines  taut.  The  samples  were  then  taped  to  the  fishing  lines,  effectively  levitating  the 
samples  in  front  of  the  detector. 

Originally,  the  samples  were  taped  to  a  block  of  polyethylene  and  were  placed  on 
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the  same  lab  jack  as  the  source.  This  approach  suffered  from  the  vibration  of  the  de¬ 
tector  cooler.  The  vibration  caused  the  block  to  skate  around  on  top  of  the  laboratory 
jack.  The  feet  of  the  girder  assembly  were  far  enough  away  from  the  detector  not  to 
be  affected  by  the  vibration,  and  the  feet  of  the  girder  assembly  are  padded.  Also,  by 
hanging  the  samples,  the  Compton  scatter  events  in  the  polyethylene,  which  would 
add  noise  to  the  results,  are  eliminated.  The  wooden  block  and  the  steel  shielding 
block  also  suffered  translation  over  long  periods  of  time  due  to  the  vibration  of  the 
cooler.  These  were  secured  to  the  jack  with  screws.  The  detector  system  has  built-in 
mechanisms  to  dampen  the  effects  of  the  cooler  vibration  on  the  detector  crystal. 
Thus,  with  the  hanging  sample  configuration,  the  sample  and  the  detector  crystal 
can  be  considered  stationary  relative  to  each  other. 


PHDs 

Det. 


Phantom 


Steel 

Shield 


Figure  27.  The  basic  setup  used  to  begin  to  characterize  the  MCST  forward  problem, 
without  collimation. 

If  Cd-109  is  being  used  for  the  source,  then  lead  cannot  be  used  for  shielding. 
Cd-109  emits  an  88  keV  7-ray.  There  is  a  good  probability  that  a  7-ray  with  88  keV 
of  energy  will  knock  out  a  K-shcll  electron  in  the  lead.  When  the  vacancy  then  is 
filled  with  another  electron,  an  X-ray  between  72  keV  and  88  keV  is  emitted  [20]. 
The  measured  backscattered  photons  from  the  sample  are  between  about  65  keV 
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and  80  keV.  Essentially  the  lead  shielding  becomes  an  additional,  unwanted  source, 
emitting  photons  in  the  energy  range  of  interest.  A  lower  Z  material  with  lower  K- 
shell  transition  energies  was  chosen  for  shielding.  Steel  was  used  in  this  case,  based 
on  availability.  A  5  cm  thick  block  was  used,  sufficiently  stopping  direct  photon  travel 
from  the  source  to  the  detector. 


Figure  28.  The  MCST  setup,  viewed  from  the  front,  looking  back  at  the  detector. 


Alternatively,  a  higher  energy  source  could  be  used,  such  that  the  photons  scat¬ 
tered  in  the  sample  would  be  easily  distinguishable  from  the  lead  X-rays.  Cd-109 
was  chosen  because  it  was  the  most  active  source  available  at  the  time  and  because  it 
was  in  the  smallest  package,  making  the  point  source  approximation  more  legitimate. 
Brian  Evans  et  al  had  originally  found  the  88  keV  photons  from  Cd-109  to  be  opti- 
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mal  for  Compton  backscatter  studies  in  metals  for  several  reasons  [25].  First,  Cd-109 
does  not  have  significant  7-ray  emissions  of  energies  higher  than  88  keV,  which  would 
downscatter  into  the  energy  range  of  interest,  introducing  noise.  Second,  the  energy 
must  be  high  enough  so  that  Compton  scatter  is  the  dominant  mode  of  interaction 
in  the  material.  Third,  the  energy  must  be  low  enough  such  that  the  mean  free  path 
in  the  material  is  not  too  high.  The  higher  the  mean  free  path,  the  less  chance  of 
interaction  there  will  be  in  the  sample. 

In  millimeters,  based  on  the  coordinate  system  shown  in  Figure  15,  the  source 
was  placed  at  (122.4  +/-  1,  175  +/-  1,  52.5  +/-  2.5).  The  sample  was  placed  at  205.4 
mm  +/-  1  mm  in  the  x  direction.  Typically  in  the  y  and  z  directions,  the  source  was 
located  at  80  and  50  mm,  respectively,  though  these  value  varied  by  centimeters  from 
test  to  test. 

Numerous  samples  were  fabricated  by  the  AFIT  model  shop.  Each  sample  was 
about  5  mm  x  5  mm  x  0.7  mm.  Samples  were  made  of  copper,  aluminum,  brass  and 
steel.  A  CuCu  “phantom”  was  constructed,  in  which  two  of  the  copper  samples  were 
arranged  5  cm  from  each  other,  center  to  center,  and  glued  to  stiff  paper.  Another 
phantom  was  made  in  which  one  sample  of  aluminum  and  two  samples  of  copper 
were  arranged  in  an  isosceles  triangle  pattern  (Figure  20).  The  copper  strips  were 
spaced  2  cm  apart  and  the  aluminum  sample  was  placed  2  cm  above  the  copper  strips, 
centered. 

4.4  Coded  Aperture  Experimental  Setup 

For  the  coded- aperture  experiments,  the  coordinate  system  shown  in  Figure  15  was 
used  again.  The  programmable  collimator  was  placed  on  a  laboratory  jack  as  close  to 
the  detector  in  the  x  direction  as  was  possible.  By  minimizing  the  detector-collimator 
distance,  the  magnification  parameter  was  also  minimized.  In  the  x  direction,  mea- 
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Figure  29.  The  basic  setup  for  the  coded-aperture  experiments.  Distances  used  are 
given  in  Table  3.  Notice  that  collimator  depth  is  significant  with  the  distances  used. 

surements  were  made  from  the  middle  of  the  detector  crystal  and  the  middle  of  the 
programmable  collimator. 

Early  on  it  was  discovered  that  the  PEEK  pings  did  not  move  easily  in  the  cham¬ 
bers,  even  when  liberal  amounts  of  way  oil  were  used  as  lubrication.  The  plugs  could 
be  pushed  in  to  the  “on”  position  easily  enough,  but  problems  occurred  when  pulling 
the  plugs  out  to  the  “off”  position.  Sometimes  excessive  force  would  be  needed  to 
pull  the  plugs  out.  Each  plug  had  been  machined  a  bit  larger  than  the  chambers, 
and  was  broached  so  as  to  form  the  tightest  possible  seal.  A  few  of  the  plugs  fit 
too  tightly.  Sometimes  in  dry  tests,  when  pulling  a  given  plug  out,  the  plug  would 
be  pulled  out  of  the  chamber  entirely.  The  experimenter  would  have  to  apply  great 
strength  to  overcome  the  static  friction.  Once  the  plug  started  moving  he  would  not 
be  able  to  stop  himself  in  time,  and  the  plug  would  come  completely  out.  If  this  were 
to  happen  during  an  experiment,  molten  AIM-70  would  gush  from  the  open  chamber 
onto  the  experimenter  and  the  setup. 

The  experimental  setup  was  adjusted  to  temporarily  circumvent  the  above  short¬ 
comings  of  the  device  until  a  permanent  solution  is  found.  To  change  the  mask 
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Table  3.  Source  Positions  for  the  Coded-Aperture  Experiment 


Mask-Det  Dist  [cm] 

Src-Det  Dist  [cm] 

Magnification 

Position  0 

5.6 

30.3 

1.33 

Position  1 

5.6 

60.6 

1.13 

Position  2 

5.6 

177.6 

1.04 

Position  F 

8 

60.6 

1.16 

configuration,  the  programmable  collimator  was  taken  to  a  fume  hood.  The  device 
was  drained  of  AIM-70  before  the  configuration  was  changed.  After  the  configuration 
change,  the  collimator  was  refilled  with  AIM-70.  The  collimator  was  then  placed 
back  in  to  the  experimental  setup.  This  method  obviously  is  not  ideal  and  takes 
much  longer  than  the  intended  method.  However,  one  advantage  of  this  method  is 
that  the  collimator  can  cool  before  being  used  in  the  experiment.  There  had  been 
concern  about  a  hot  programmable  collimator  being  too  close  to  the  detector.  The 
detector  crystal  must  be  kept  at  about  50  K.  By  letting  the  AIM-70  cool  and  solidify, 
one  need  not  worry  about  heating  up  the  detector.  Also,  the  additional  complexities 
of  having  a  hot  plate  in  the  experimental  area,  along  with  the  need  to  monitor  the 
molten  AIM-70,  are  eliminated. 


Figure  30.  Three  selected  sources  used  in  the  experiments.  Left:  A  115  /iCi,  Co-57 
button  source.  Center:  An  intermittent  38.13  /zCi,  Co-57  Rod.  Right:  A  47.66  /jCi, 
Co-57,  Flexible  Rope 
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Three  source  positions  were  used  in  order  to  try  different  magnifications.  Various 
sources  were  imaged.  The  Co-57  sources  were  used  most  in  the  coded-aperture  testing. 
The  button  source  was  useful  for  approximating  a  point  source.  The  Co-57  rod  is 
about  15  cm  long,  and  is  “dotted,”  or  intermittent.  There  are  1  cm  long  lead  pieces 
spaced  1  cm  apart  down  the  length  of  the  rod.  The  Co-57  rope  was  arranged  in 
a  distinctive  shape.  For  each  mask  configuration,  the  sources  were  placed  at  the 
different  positions  and  data  were  taken  for  various  amounts  of  time. 

Because  the  detector  is  energy-sensitive,  an  energy  window  could  be  set.  For 
instance,  with  a  Co-57  source,  the  energy  window  is  set  such  that  the  only  photon 
events  recorded  in  the  image  are  those  with  energies  of  122.1  keV  +/-  3  keV  or 
136.5  keV  +/-  3  keV.  By  limiting  the  energy  of  the  photons  allowed,  photons  that 
have  Compton  scattered  can  be  excluded.  Indeed,  only  photons  that  pass  straight 
through  a  chamber  in  the  “on”  configuration  are  desired.  The  energy  sensitivity  of 
the  detector  helps  to  exclude  those  photons  that  scatter  in  closed  chambers.  The 
energy  sensitivity  also  serves  to  limit  the  background  y-rays  that  are  measured,  thus 
decreasing  the  noise  in  the  raw  data.  Even  so,  the  detector  was  surrounded  on  the 
back  and  sides  with  walls  of  lead  bricks,  to  shield  the  detector  from  scattered  photons 
originating  from  other  experiments  taking  place  in  the  laboratory. 


Table  4.  A  selection  of  the  sources  used. 


Designator 

Isotope 

Form 

Activity  [/iCi] 

Ref.  Date 

7  Energy  [keV] 

T-122 

Co-57 

Button 

115 

15NOV11 

122.06,  136.47 

Rod 

Co-57 

Rod 

38.13 

04JAN12 

122.06,  136.47 

Rope 

Co-57 

Rope 

47.66 

04JAN12 

122.06,  136.47 

T-124 

Co-60 

Planchet 

5.85 

04JAN12 

1173,  1333 

T-098 

Cd-109 

Button 

101.8 

07NOV11 

88.03 
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V.  Results  and  Discussion 


5.1  Coded- Aperture  Image  Reconstruction  Using  Decorrelation  Methods 

A  good  first  place  for  data  analysis  is  with  the  correlation  methods.  Historically 
these  methods  were  used  first  and  are  still  used  today  in  some  applications.  They 
are  much  less  computationally  intensive  and  much  easier  to  implement.  As  for  the 
choice  of  the  matched  filter,  there  are  a  few  different  options.  One  is  to  use  the  mask 
pattern  as  the  matched  filter  G.  Another  choice  is  the  Brown  matched  filter,  where 
the  zeros  are  replaced  with  “-l’s.” 


The  Matched  Filter 


2  4  6  8  10 


The  Raw  Data  from  the  Detector,  Colorbar  in  Counts 
111116  Co57Point  Posl  1200sec 


50  100 

Subpixel  Column 


Object  Estimate,  Decorrelated  w/  Same  Filter  Object  Estimate,  Decorrelated  w/  BrownFilter  Object  Estimate,  Decorrelated  w /  BrownFilter,  Negative  Removed 

Colorbar  is  Probability  Density  lmA-21  Colorbar  is  Probability  Density  Im^l  Colorbar  is  Probability  Density  [mA-2] 


Figure  31.  Images  from  a  test  using  the  Co-57  button  source  at  position  1,  over  1200 
s.  Top-Left:  The  mask  pattern.  Top-Right:  The  raw  data.  Bottom-Left:  The  object 
estimate,  using  the  mask  pattern  as  G.  Bottom-Center:  The  object  estimate  using 
a  Brown  matched  filter.  Bottom-Right:  The  object  estimate  using  a  Brown  matched 
filter,  setting  negatives  to  zero. 
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A  Matlab®  script  was  written  to  decorrelate  the  raw  data.  The  mask  pattern 
is  input  as  a  matrix,  as  well  as  the  directory  of  the  raw  data  text  hie.  The  routine 
imports  the  data,  storing  it  as  a  two-dimensional  matrix.  Next,  the  routine  over¬ 
resolves  both  the  raw  data  matrix  and  the  mask  matrix  by  the  same  factor  so  that 
they  maintain  equal  scaling.  The  pattern  size  is  increased  by  a  factor  of  m,  where 
m  is  the  magnification  expected.  Because  the  matrices  are  over-resolved,  the  mask 
shadow  can  be  rounded  to  a  much  finer  range  of  sizes.  By  including  the  effects  of 
magnification,  the  scale  of  the  decorrclated  image,  or  object  estimate,  is  correct. 

After  the  over-resolved  raw  data  and  G  matrices  are  defined,  the  script  invokes 
the  built-in,  cross-correlation  function,  such  that 

d  =  P*G  (21) 

where  P  is  the  raw  data  and  O  is  the  object  estimate.  The  numerical  cross  correlation 
function  outputs  a  matrix  that  is  bigger  than  the  two  matrices  that  are  input,  as 
expected.  If  P  is  of  size  \Mp,Np]  and  G  is  of  size  [MG,NG],  then  O  will  be  of  size 
[. MP  +  Mg  —  1,  NP  +  Ng  —  1],  Abstractly,  the  cross  correlation  process  can  be  thought 
of  as  discretely  sliding  G  over  P  and  recording  how  well  the  two  “agree”  at  each 
location.  In  this  case  only  the  center  part  of  O  was  used.  More  specifically,  the  center 
[MP,Np\  square  of  O  was  extracted,  and  the  rest  was  ignored.  The  outside  portions 
do  have  meaning,  representing  instances  where  the  mask  shadow  is  not  fully  on  the 
detector  crystal.  However,  it  was  found  that  the  outside  portions  were  not  reliable, 
and  so  they  were  discarded.  Furthermore,  the  sources  were  placed  in  a  narrow  enough 
solid  angle  so  that  the  shadow  would  not  project  off  of  the  side  of  the  detector  crystal. 
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The  script  also  normalizes  O.  That  is, 


O'  = 


(22) 


By  normalizing  the  object  estimate,  the  units  of  O'  become  cm-2.  The  values  of  O' 
are  then  values  of  probability  density.  If  O'  were  integrated  over  its  limits  the  result 
would  be  unity.  The  source  is  assumed  to  be  located  entirely  somewhere  within  the 
spatial  limits  of  the  O'  matrix  and  O'  provides  the  probability  of  finding  it  at  a  given 
location  within  those  limits. 

Regarding  inherent  noise  and  filtering  artifacts,  perhaps  the  most  instructive  of 
the  tests  performed  were  those  that  used  the  Co-57  button  source.  In  the  bottom-left 
portion  of  Figure  31,  one  can  see  the  results  of  using  d  as  G.  A  large  amount  of 
inherent  noise,  or  artifacts  that  arise  from  the  decorrelation  process,  are  seen.  This 
inherent  noise  is  manifested  throughout  the  reconstructed  image.  Though  the  button 
source  is  still  clear  in  this  case,  it  may  not  be  when  short  detection  times  are  used  or 
when  a  weak  source  is  imaged.  The  Brown  matched  filter  is  able  to  smooth  out  the 
inherent  noise,  to  a  somewhat  constant  baseline,  as  seen  in  the  bottom-center  image. 
Nonuniformities  are  certainly  present,  but  the  source  is  more  easily  distinguishable. 

A  peculiarity  that  arises  when  using  the  Brown  type  matched  filter  is  that  O' 
can  take  on  negative  values.  Obviously  a  negative  probability  is  not  meaningful 
physically.  One  way  to  account  for  these  negative  values  is  to  throw  them  away 
entirely.  Assuming  the  negative  values  of  probability  to  be  non-physical,  the  locations 
corresponding  to  those  negative  values  are  given  zero  probability  of  containing  the 
source.  This  explanation  may  not  sit  well  with  some.  Another  way  to  think  about 
this  method  is  that  the  baseline  is  being  raised.  If  the  Brown  matched  filter  causes 
a  relatively  constant  level  of  inherent  noise  across  0,  this  constant  noise  simply  can 
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be  removed.  Some  signal  may  be  lost  but  the  SNR  should  increase.  Viewing  the 
bottom- right  image  in  Figure  31,  the  advantage  of  this  correction  is  clear.  In  practice 
one  could  move  the  threshold  or  baseline  even  higher,  discarding  even  more  of  the 
noise. 

It  is  known  that  coded-apertures  are  not  as  effective  for  extended  source  imaging. 
This  is  why  astronomy  became  the  first  real  use  for  the  coded-aperture.  Stars  are 
essentially  point  sources.  Nevertheless,  the  rod  and  the  rope  sources  were  used  to 
characterize  the  programmable  collimator’s  ability  to  image  extended  sources.  One 
can  see  a  bit  of  the  rod  in  Figure  32.  However,  the  noise  is  much  more  pronounced 
in  this  case.  In  Figure  33  the  rope  is  mostly  indistinguishable. 

There  are  a  few  factors  that  worsen  the  object  estimate  in  this  particular  case. 
First,  the  detector  crystal  is  round,  as  can  be  seen  in  the  left  image  of  Figure  32. 
Not  only  is  data  lost  because  of  this  fact,  but  data  is  lost  non-uniformly  in  space. 
The  decorrelation  methods  do  not  have  a  good  way  to  account  for  this  fact,  although 
ML-EM  does. 

Secondly,  the  total  projection  through  the  collimator  can  become  bigger  than  the 
detector  crystal  if  there  is  any  significant  magnification.  Again,  this  causes  a  spatially 
dependent  loss  of  data  that  the  decorrelation  methods  are  ill-suited  to  handle.  One 
possible  solution  is  to  have  either  a  bigger  detector  or  a  smaller  collimator  so  that 
the  projection  stays  on  the  crystal,  away  from  the  corners.  A  bigger  detector  crystal 
is  not  foreseeable  with  current  equipment  and  funds.  A  smaller  collimator  could  be 
fabricated  given  two  to  five  months.  Another  solution  is  to  move  the  source  farther 
away  from  the  collimator,  decreasing  the  magnification.  This  has  the  disadvantage 
of  lower  count  rates.  Also,  given  the  resulting  smaller  solid  angle  subtended  by  the 
source,  from  the  detector,  the  resolution  of  the  object  image  would  be  worse  if  the 
source  were  farther  away. 
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The  Raw  Data  from  the  Detector,  Colorbar  in  Counts  Object  Estimate,  Decorrelated  w/  BrownFilter  Object  Estimate,  Decorrelated  w/  BrownFilter,  Negatives  Removed 


Subpixel  Column  [cm]  [cm] 


Figure  32.  Data  from  the  rod  at  position  zero.  Left:  The  raw  data.  Center:  The  object 
estimate,  using  a  Brown  matched  filter.  Right:  The  object  estimate,  Brown  matched 
filter,  negative  values  set  to  zero. 


Thirdly,  vignetting  was  significant  with  the  dimensions  used  in  this  study.  As 
mentioned  in  Section  3.1,  there  is  an  analytic  formula  which  can  correct  the  vignetting 
effects  [50].  Such  a  correctional  factor  could  be  applied  in  the  future.  Essentially,  the 
raw  data  pixels  would  be  weighted  based  on  their  location,  where  the  center  pixels 
would  have  a  higher  weight  than  the  edge  pixels. 


Co-57  Rope,  OHat,  Brown  Filter,  Co-57  Rope,  OHat,  Brown  Filter,  No  Negs 

Colorbar  is  Probability  Density  [cmA-2]  ¥  in'3  Colorbar  is  Probability  Density  [cmA-2]  x10~3 


[cm]  [cm] 

Figure  33.  Left:  The  Co-57  rope  source.  Center:  The  object  estimate  using  a  Brown 
matched  filter.  Right:  The  object  estimate,  Brown  matched  filter,  negative  values  set 
to  zero. 
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5.2  Coded- Aperture  Image  Reconstruction  Using  ML-EM 

The  same  data  analyzed  with  correlation  techniques  in  Section  5.1  were  also  ana¬ 
lyzed  using  the  ML-EM  method.  The  matrix  A  was  found  as  described  in  Section  3.1. 
A  script  was  written  to  import  the  data  from  the  detector  and  to  run  the  algorithm. 
The  ML-EM  equation  for  imaging,  reproduced  here  from  Section  2.4,  is  [53,  40,  62]: 


Afd 

^new  _  3 


y  a 


Vi 


(23) 


Ei  A  ~  Efc  AfcA*;  +  b 
The  method  performed  well  with  point  sources.  The  number  of  iterations  required 
for  a  satisfactory  object  image  varies.  When  imaging  a  point  source,  the  source 
became  distinct  with  just  a  few  iterations.  Figure  34  shows  the  object  estimate  of 
a  Co-57  point  source.  After  just  4  iterations,  the  source  is  distinct,  though  with  a 
halo.  After  200  iterations,  the  source  is  reduced  to  a  circle  with  a  diameter  of  about 
1  cm.  The  button  source  has  a  physical  diameter  of  1.5  cm.  After  1400  iterations, 
the  pixels  of  the  image  become  more  differentiated.  The  right  image  of  Figure  34  is 
an  example  of  ML-EM  convergence.  After  14,000  iterations,  for  instance,  the  image 
would  look  similar  to  that  of  1,400  iterations. 


ML-EM  Src  Plane  Est.  Iterations:  4  ML-EM  Src  Plane  Est.  Iterations:  200  ML-EM  Src  Plane  Est  Iterations:  1400 


Figure  34.  The  Co-57  button  source  was  imaged  for  30  seconds.  Shown  is  the  object 
estimate  after  4  iterations  (Left),  200  iterations  (Center)  and  1400  iterations  (Right). 
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The  average  amount  that  each  image  pixel  changed  from  iteration  to  iteration  is 
shown  in  Figure  35.  Notice  the  drastic  changes  below  about  100  iterations.  Indeed, 
convergence  can  be  assumed  above  a  few  hundred  iterations,  where  the  plots  become 
asymptotic.  One  should  note  that  full  convergence  in  ML-EM  often  results  in  some 
non-physical  artifacts.  These  non-physical  artifacts  arise  due  to  limitations  in  the 
forward  model  predications  as  well  as  from  limited  sampling  in  the  raw  data.  In  the 
right-hand  image  of  Figure  34,  for  instance,  the  blotchy  pixelation  shown  is  probably 
not  physical.  In  a  commercial  system  the  number  of  iterations  would  be  limited,  in 
order  to  prevent  these  artifacts  from  developing.  By  limiting  the  number  of  iterations, 
the  reconstructed  image  may  not  reach  the  same  spatial  resolution,  though  the  number 
of  non-physical  artifacts  is  less.  Both  fully  converged  and  partially  converged  images 
are  presented  in  this  document. 
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Figure  35. 
next.  Left: 
rod. 


The  average  change  in  reconstructed  image  pixels  from  one  iteration  to  the 
Data  for  the  Co-57  button  source.  Right:  Data  for  the  Co-57  intermittent 


The  field-of-view  was  constrained  because  of  the  collimator’s  depth.  It  was  dis¬ 
covered  that  the  ML-EM  method  would  fail  if  the  discretized  source  plane  did  not 
fit  within  the  field-of-view.  Specifically,  the  reconstructed  image  would  be  populated 
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entirely  with  NaNs.  This  outcome  is  related  to  the  collimator’s  depth.  As  described 
in  Section  3.1,  the  forward  model  finds  a  vector  between  each  source  pixel  and  all 
of  the  detector  pixels.  If  a  photon  traveling  along  a  vector  from  source  pixel  j  to 
detector  pixel  i  is  blocked  by  the  collimator,  then  Al3  is  zero.  Now  suppose,  as  an 
example,  that  source  pixel  j  =  2  is  completely  outside  of  the  collimator’s  field-of- 
view.  In  this  case,  photons  traveling  from  source  pixel  2  to  any  detector  pixel  are 
blocked.  Therefore  the  term  JA  A ^  goes  to  zero.  This  term  is  located  in  the  de¬ 
nominator  of  Equation  23.  In  computers,  division  by  zero  results  in  an  assignment  of 
the  value  NaN.  Continuing  with  the  example,  A^6"  thus  would  go  to  NaN.  The  NaN 
propagates  to  other  j  positions  in  A  because  of  the  ^fcAjfcAfc  term.  NaN  times  a 
number  equals  NaN.  Within  a  couple  of  iterations,  A  is  filled  completely  with  NaNs. 
Source  configurations  must  be  chosen  in  which  the  source  is  completely  within  the 
ficld-of-view. 

The  parameter  b  represents  the  average  number  of  counts,  per  pixel,  that  are  not 
attributable  to  unscattered  rays  from  the  source.  That  is,  b  represents  the  average 
number  of  counts  that  are  not  of  the  desired  signal.  If  this  parameter  is  set  too 
high,  true  data  can  be  lost.  Less  active  regions  of  the  source  plane  are  set  to  zero 
artificially.  On  the  other  hand,  if  b  is  set  too  low,  noise  and  background  counts  corrupt 
the  object  estimate.  The  number  of  iterations  and  b  were  set  qualitatively,  based  on 
how  the  object  estimate  looked  after  the  script  was  finished.  Typically  the  number  of 
iterations  needed  for  convergence  ranged  between  300  and  800.  b  was  usually  set  to 
be  10%  of  the  counts  in  the  raw  data  pixel  with  the  most  counts.  A  more  quantitative 
means  for  determining  these  two  values  would  be  useful  in  a  commercialized  system. 

The  data  taken  using  the  extended  sources  were  analyzed  using  ML-EM.  The 
method  performed  well,  especially  as  compared  with  the  correlation  methods.  As 
seen  in  Figure  36,  the  full  structure  of  both  the  rod  and  the  rope  is  visible.  The 
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images  do  not  hone  in  exactly  to  the  width  of  the  rod  or  the  rope.  The  radioactive 
material  in  both  the  rod  and  the  rope  has  a  width  of  2  mm.  In  the  reconstructed 
images  shown  in  Figure  36,  the  rod  has  a  width  ranging  from  3-6  mm  and  the  rope 
has  a  width  ranging  from  4-5  mm.  This  limited  resolution  does  not  stem  from  any 
problem  with  the  ML-EM  method  but  rather  is  due  to  an  interplay  between  the 
magnification,  the  cross-sectional  area  of  a  given  collimator  chamber,  and  the  area  of 
a  given  detector  pixel.  The  magnification  of  the  system  could  be  altered  by  changing 
the  source-mask  and  mask-detector  distances.  If  the  magnification  is  increased,  the 
resolution  of  the  reconstructed  image  increases  but  then  less  of  the  source  plane  lies 
within  the  system’s  field-of-view,  and  vice-versa.  The  magnification  and  the  resolution 
are  directly  related.  The  magnification  was  chosen  in  this  case  so  as  to  fit  a  significant 
portion  of  the  extended  sources  within  the  field-of-view. 


Co-57  Rod,  ML-EM  Src  Est.  Iterations:  800  Co-57  Rope,  ML-EM  Src  Est.  Iterations:  800 

Colorbar  is  Probability  Density  [cm ^2]  Colorbar  is  Probability  Density  [cm ^2] 


Figure  36.  Top-Left:  The  Co-57  intermittent  rod.  Bottom-Left:  ML-EM  reconstruction 
of  the  rod  after  800  iterations.  Data  taken  for  1  hr.  Top-Right:  The  Co-57  rope. 
Bottom-Right:  ML-EM  reconstruction  of  the  rope  after  800  iterations.  Data  taken  for 
1  hr. 
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The  primary  advantage  of  the  ML-EM  method  over  correlation  methods  is  its 
ability  to  incorporate  certain  effects  and  physical  characteristics  of  the  detector  that 
cannot  be  accounted  for  using  correlation  methods.  For  instance,  the  A  matrix  takes 
into  account  the  fact  that  every  tenth  subpixel  in  the  detector  is  dead,  and  that 
the  detector  is  round,  not  square.  And,  as  mentioned  before,  the  forward  problem 
results  naturally  compensate  for  vignetting  effects.  With  the  continual  improvement 
of  computer  memory  and  speed,  ML-EM  will  only  become  better  as  compared  to 
the  correlation  methods.  For  this  problem,  algorithm  computation  times  were  trivial 
except  when  the  number  of  source  plane  pixels,  M,  in  the  matrix  B  were  greater  than 
about  20,000,  in  which  case  computation  times  were  on  the  order  of  ten  minutes.  A 
personal  computer  was  used. 

5.3  Quantitative  Image  Analysis  and  SNR 

Comparing  two  images,  the  human  eye  is  usually  good  at  determining  when  an 
image  is  better  than  another.  In  an  image  of  a  radioactive  source,  the  source  should  be 
bright  (high  values)  and  sharply  defined  whereas  the  rest  of  the  image  should  be  black 
(values  of  zero)  and  free  of  any  artifacts.  The  eye  can  pick  these  sorts  of  things  out,  but 
a  more  quantitative  analysis  is  almost  always  preferable  in  science  and  engineering. 
Two  methods  for  finding  the  Signal-to- Noise  Ratio  (SNR)  in  reconstructed  images 
were  investigated. 

Both  methods  used  the  definition  of  SNR  from  Cunningham  et  al  as  a  starting 
point  [19].  The  size  and  location  of  the  source  being  imaged  is  known.  The  source, 
being  the  signal,  is  removed  from  the  reconstructed  image.  Only  noise  remains  in  the 
image.  The  signal  and  the  noise  then  can  be  compared. 

A  script  was  written  to  determine  the  SNR  in  images  of  button  sources.  The  user 
is  able  to  define  the  center  and  the  radius  of  a  circle  which  represents  the  button 
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Figure  37.  From  the  object  estimate  (left),  a  circle  of  signal  is  found,  by  hand,  leaving 
the  rest  to  be  noise  (right). 


source.  The  values  of  the  pixels  in  the  image  that  are  covered  by  the  circle  are  placed 
in  a  length  Mb  vector  called  B,  where  Mb  is  the  number  of  pixels  that  are  covered  by 
the  user-defined  circle.  The  values  of  the  pixels  not  covered  by  the  circle  are  placed 
into  a  vector  called  C  of  length  Me ,  where  Me  is  the  number  of  pixels  that  are  not 
covered  by  the  circle.  Notice  that  Mb  +  Me  equals  the  number  of  pixels  in  the  image. 

In  the  first  version  of  SNR  used,  the  mean  of  the  values  in  the  signal  list,  B ,  and 
the  mean  of  the  values  in  the  noise  list,  C,  are  found.  The  SNR  is  simply: 

SNR  =  =  (24) 

c 

This  formulation  works  well,  except  in  the  correlation  method  that  uses  the  Brown 
matched  filter.  When  the  Brown  matched  filter  is  used,  the  object  estimate  contains 
negative  values.  The  negative  noise  is  just  as  detrimental  to  the  image  quality  as  is 
the  positive  noise.  However,  if  a  simple  mean  is  taken  of  the  noise  list,  then  many  of 
the  negative  noise  values  would  cancel  out  the  positive  values.  C  would  be  close  to 
zero,  yielding  an  artificially  high  SNR.  Another  version  of  the  SNR  compensates  for 
these  negative  values  by  using  the  Root  Mean  Square  (RMS)  of  the  signal  list  and  of 
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the  noise  list.  The  RMS  signal,  S,  is  defined  as 


S 


(25) 


and  the  RMS  noise,  £  is 

I  sr^Mc  ~7p 

C  =  V  Mb  (26) 

The  SNRrms  is  defined  as 

SNRrms  =  |  (27) 

The  signal  and  noise  distributions  were  examined  for  general  insight  and  to  de¬ 
termine  if  an  appropriate  simplification  could  be  made.  An  example  is  shown  in 
Figure  38.  The  noise  that  is  found  in  images  reconstructed  with  correlation  methods 
generally  adheres  to  the  Poisson  distribution,  even  through  the  decoding  process.  Ra¬ 
dioactive  counting  experiments  in  which  the  counting  time  is  short  compared  to  the 
half  lives  of  the  nuclides  used  will  follow  Poisson  statistics.  In  this  case,  the  standard 
deviation  is  equal  to  the  square  root  of  the  mean  of  the  distribution  [37,  64], 
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Figure  38.  Left:  A  histogram  of  C  taken  from  an  image  of  the  Cd-109  source  with 
background.  Right:  A  histogram  of  B  from  an  image  of  the  Cd-109  source  with  back¬ 
ground. 
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The  signal  list,  B,  from  the  correlation  method  images,  as  well  as  both  B  and  C 
from  the  ML- EM  images,  do  not  adhere  to  any  known  distribution.  For  thoroughness, 
the  standard  deviation  for  B  and  C  was  calculated  outright  for  all  images.  For  the 
first  SNR  formulation, 


Under  the  SNRrms  formulation, 

(30) 


In  the  calculation  of  SNR,  the  variances  add  [68].  Eventually  the  images  could 
be  conditioned,  removing  noise  without  significant  loss  of  the  true  signal.  The  goal 
of  this  research  is  not  to  produce  the  best  possible  images  but  to  show  the  relative 
improvement  between  images  when  multiple  mask  patterns  are  used. 

5.4  Image  Improvement  from  a  Multiple  Mask  Sequence 

As  a  coded-aperture,  the  overarching  goal  of  the  programmable  collimator  is  better 
image  quality.  This  is  accomplished  by  using  different  mask  patterns,  which  provide 
semi-independent  views  of  the  source.  The  different  views  provide  new  data  that 
helps  to  cancel  out  noise  and  irregularities  resulting  from  the  previous  views’  data. 
The  programmable  collimator  can  improve  imaging  systems  in  two  different  fashions. 
It  can  either  make  better  images  in  the  same  measurement  time  as  compared  to  a 
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traditional  coded-aperture,  or  it  can  arrive  at  the  same  quality  image  that  a  traditional 
coded-aperture  would  produce  but  in  less  time.  The  former  situation  was  explored 
in  this  test  campaign.  The  Co-57  button  source,  the  Cd-109  button  source  (with 
background  injection),  the  itermittent  rod,  and  the  rope  were  all  imaged  at  Position 
F  (See  Tables  4  and  3).  The  rope  was  only  imaged  with  up  to  five  masks,  rather  than 
ten. 
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Figure  39.  The  10-mask  sequence  of  random  masks  used  for  testing.  The  sequence  was 
found  by  the  program  described  in  Section  3.3. 


The  Cd-109,  rod,  and  rope  sources  were  all  imaged  for  a  total  of  1  hour.  The  Co-57 
source  was  imaged  for  a  total  of  5  minutes.  As  more  masks  were  used,  the  total  time 
was  divided  between  each  mask.  For  example,  when  the  intermittent  rod  was  viewed 
through  two  masks,  data  were  taken  with  each  mask  for  30  minutes.  When  three 
masks  were  used,  data  were  taken  through  each  mask  for  20  minutes,  and  so  forth. 
In  this  way  the  total  measurement  time  remained  the  same  so  that  fair  comparisons 
could  be  made  between  sequences  with  different  quantities  of  masks. 

To  reduce  the  number  of  mask  changes  to  be  made  to  the  collimator,  all  of  the 
sources  were  viewed  through  each  mask  before  changing  the  mask  pattern.  All  re¬ 
quired  data  collection  times  were  found  for  each  mask  before  changing  the  pattern. 
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Unfortunately,  the  proprietary  HPGe  detector  software  did  not  allow  automation, 
meaning  that  each  data-run  was  initiated  manually.  If  automation  had  been  an  op¬ 
tion,  more  more  mask  patterns  would  have  been  used.  Nevertheless,  good  results 
were  found  with  ten  mask  patterns. 

A  standard  coded-aperture  setup  was  used  for  the  test  campaign  (Figure  29).  The 
position  of  the  collimator  and  the  sources  were  marked  on  their  respective  stands 
using  a  sharp  pen.  The  collimator  was  taken  to  a  different  laboratory  with  a  fume 
hood  for  reconfiguration.  When  returned,  it  was  positioned  according  to  the  marks 
on  the  stand.  The  sources  were  repositioned  similarly  between  tests.  This  method 
led  to  some  net  movement  or  magnification/demagnification  of  the  source  in  the 
reconstructed  images,  as  explained  in  Subsection  5.4.1.  Images  were  reconstructed 
using  the  Brown  matched  filter  decorrelation  technique  and  using  ML-EM.  When 
multiple  masks  were  used,  a  reconstruction  technique  was  performed  on  the  data  from 
each  mask  pattern.  The  reconstructed  images  were  weighted  equally  and  summed 
together.  The  summed  image  was  then  normalized  to  allow  for  legitimate  comparisons 
between  different  sequences.  Image  improvement  was  demonstrated  with  all  sources, 
as  described  in  the  following  subsections. 


Table  5.  The  data-runs  performed  for  the  Cd-109  source.  Each  “1”  represents  a  single 
data-run. 


Mask 

3600  s 

1800  s 

1200  s 

900  s 

720  s 

600  s 

514  s 

450  s 

400  s 

360  s 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

2 

1 

1 

1 

1 

1 

1 

1 

1 

1 

3 

1 

1 

1 

1 

1 

1 

1 

1 

4 

1 

1 

1 

1 

1 

1 

1 
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1 

1 

1 

1 

1 

1 
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1 

1 

1 

1 
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1 

1 

1 

1 
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1 

1 

1 

9 

1 

1 

10 

1 
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5.4.1  Image  Shift  Correction 


It  was  discovered  that  the  collimator  had  translated  between  mask  patterns.  When 
the  mask  was  reconhgnred  in  a  different  laboratory,  it  was  not  returned  to  the  exact 
same  position  in  the  imaging  laboratory.  This  caused  a  shift  in  the  reconstructed 
images. 


Table  6.  The  coordinates  of  the  center /maximum  pixel  of  the  reconstructed  button 
source  for  each  of  the  ten  masks,  in  mm.  Co-57  dataset  1. 


Mask  Num. 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

Row 

48 

50 

50 

51 

50 

49 

50 

50 

52 

51 

Column 

44 

48 

51 

50 

46 

52 

46 

51 

51 

51 

For  every  mask  pattern,  data  were  taken  using  the  Co-57  button  source.  The  ML- 
EM  images  of  the  Co-57  source  were  so  clean  that  the  shift  could  be  discerned.  In  all 
reconstructed  images,  the  pixel  with  the  highest  probability  density  value  also  was 
the  pixel  at  the  center  of  the  reconstructed  button  source.  The  coordinates  of  these 
maximum/center  pixels  were  recorded.  The  results  from  one  trial  are  given  in  Table 
6.  Notice  that  the  center  pixel  varies  by  as  much  as  8  mm  from  mask  to  mask.  Most 
of  the  image  movement  is  horizontal,  meaning  that  the  biggest  changes  in  collimator 
placement  were  lateral. 

Results  from  multiple  trials  were  averaged,  giving  the  expected  shift  from  mask 
to  mask.  The  values  for  the  first  mask  were  used  as  the  true  position.  Images  formed 
using  subsequent  mask  patterns  were  shifted  based  on  the  difference  between  that 
mask’s  reference  coordinates  and  the  first  mask’s  reference  coordinates.  In  this  way 
the  effects  of  collimator  movement  were  reduced.  Changes  in  magnification  due  to 
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collimator  movement  were  not  accounted  for.  However,  translation  of  the  collimator 
by  a  few  millimeters  would  change  the  magnification  only  negligibly. 

For  future  experiments,  an  assembly  should  be  built  to  ensure  that  the  collimator 
does  not  change  position.  Of  course,  this  system  is  only  a  proof  of  concept.  A  real 
system  would  be  automatically  reconhgurable  using  computer  control,  eliminating 
the  need  to  ever  move  the  collimator. 

5.4.2  Co-57 

The  most  straightforward  test  was  that  which  was  performed  on  the  Co-57  button 
source.  This  active  and  small  source  served  as  a  good  approximation  of  a  point  source. 
Two  datasets  were  recorded.  About  47,000  events  inside  the  energy  window  were 
recorded  over  a  given,  total  datarun.  The  correlation  reconstruction  techniques  work 
well  when  a  point  source  is  imaged.  In  the  top-left  image  of  Figure  40,  the  source  is 
clearly  visible  when  just  one  mask  is  used,  though  significant  inherent  noise  dominates 
the  image.  When  five-mask  and  ten-mask  patterns  were  used,  the  inherent  noise  in 
the  decorrelated  image  was  reduced  significantly.  The  source  and  its  halo  stay  the 
same  size  regardless  of  the  number  of  masks  that  were  used. 

The  ML-EM  technique  produced  a  mostly  noiseless  image  with  a  single  mask 
pattern.  The  detail  in  the  source  itself  increases  when  more  mask  patterns  are  used. 
In  the  fully  converged  one-mask  image,  the  source  has  a  noticeable,  nonphysical, 
horizontal  discontinuity  and  an  artifact  slightly  above  it.  Both  features  are  not  present 
in  the  fully  converged  five-mask  image.  The  ten-mask  image  shows  a  tightened  source 
and  a  decrease  in  halo  size  versus  the  five-mask  image. 

The  SNR  was  found  for  each  mask  sequence  and  reconstruction  method  as  ex¬ 
plained  in  Section  5.3.  Results  for  the  correlation  method  are  shown  in  Figure  41. 
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Figure  40.  Images  of  the  Co-57  button  source  using  Brown  correlation,  ML-EM  with 
100  iterations,  and  ML-EM  with  500  iterations.  1,  5,  and  10  mask  sequences  are  shown. 
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SNRrms  clearly  increases  when  more  mask  patterns  are  used.  The  SNRrms  was  4.30 
+/-  0.0003  times  higher  when  ten  masks  were  used  versus  when  one  mask  was  used. 


Figure  41.  Plot  of  the  SNRrms  versus  the  number  of  masks  used.  Correlation  method. 
Co-57  button  source.  1  hour  total  measurement  time.  Error  bars  are  shown. 

One  interesting  feature  in  Figure  41  is  the  discontinuous  distribution  of  the  data 
points.  The  SNRrms  jumps  up  substantially  between  2-3  mask  sequences  and  between 
5-6  mask  sequences.  This  effect  is  explained  in  Section  5.5. 

In  the  ML-EM  images,  the  SNR  and  SNRrms  increased  by  as  high  as  nearly 
10,000.  It’s  unclear  why  the  SNR  peaks  for  the  7  and  8  mask  sequences.  The  radius 
of  the  signal  circle  was  increased,  and  the  plots  looked  similar,  meaning  that  the 
trends  are  not  caused  by  the  chosen  position  of  the  signal  circle.  SNR  analyses  on 
the  ML-EM  reconstructed  images  for  this  source  are  of  some  utility  but  should  not 
be  dwelt  upon.  The  results  were  included  for  completeness.  Already  with  a  single 
mask,  the  noise  in  the  reconstructed  image  is  very  low.  Essentially  the  black  areas  of 
the  bottom  right  image  in  Figure  40  are  zero.  The  actual  values  in  those  pixels  are 
something  very  low  (about  5xl0-43).  When  those  values  decrease,  on  average,  by  a 
factor  of  ten,  the  SNR  increases  by  a  factor  of  ten,  but  the  noise  in  those  areas  is  still 
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essentially  negligible.  Nevertheless,  the  upward  trend  in  the  data  presented  in  Figure 
42  is  undeniable. 


Number  of  Masks  Number  of  Masks 

Figure  42.  Left:  Plot  of  the  SNR  versus  the  number  of  masks  used  in  the  sequence, 
ML-EM  method,  Co-57  button  source,  1  hr  total  measurement  time.  Right:  Plot  of 
SNRrms  versus  number  of  masks  used.  Error  bars  are  shown  in  both. 


5.4.3  Cd-109  Button  Source  with  Background  Injection 

As  was  mentioned  in  Subsection  5.4.2,  the  ML-EM  reconstructed  image  of  the 
Co-57  button  source  was  already  good  with  a  single  mask  pattern.  In  order  to  better 
gauge  the  image  improvement,  a  higher  background  and  more  external  sources  of 
noise  were  desired.  In  order  to  accomplish  this,  the  Cd-109  source  was  placed  2  cm 
above  the  Co-57  button  source  (Table  4).  The  Cd-109  source  has  an  activity  of  about 
10%  of  the  Co-57  button  source,  for  the  y-rays  of  interest. 

The  unscattered  y-rays  from  the  Cd-109  source  are  observed  directly.  Some  y-rays 
from  the  Co-57  source  scatter  down  in  energy  to  the  energy  window  and  are  detected, 
also.  These  extra  photons  corrupt  the  image  of  the  Cd-109  source.  Figure  43  shows 
how  the  choice  of  energy  window  affects  the  reconstructed  image.  With  a  large  energy 
window,  more  of  the  scattered  photons  from  the  Co-57  source  are  observed,  and  vice 
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Figure  43.  The  effect  of  different  energy  windows.  The  colorscale  is  fixed.  Left:  84-92 
keV  window.  Center:  86-90  keV  window.  Right:  87-89  keV  window. 

versa.  With  a  smaller  energy  window,  more  of  the  background  is  blocked  out,  but 
the  count  rate  decreases.  Based  on  the  application,  a  choice  must  be  made  between  a 
more  stringent  energy  filter  or  higher  count  rates.  For  these  tests  an  energy  window  of 
86  keV  -  90  keV  was  used.  About  70,000  events  were  recorded  for  each  total  datarun. 

Note  that  the  Co-57  source,  seen  below  the  Cd-109  source  in  the  images  (Figure 
44),  is  somewhat  visible.  Even  though  the  scattered  photons  are  observed,  a  blurred 
version  of  the  Co-57  source  is  still  seen.  Based  on  the  lead  shielding  around  the  detec¬ 
tor,  and  the  AIM-70  within  the  housing  of  the  collimator,  these  photons  presumably 
came  through  the  open  chambers  of  the  collimator.  Compton  scattering  of  a  photon 
from  122  keV  to  88  keV  results  in  a  directional  change  of  about  130  degrees,  ruling 
out  single-scattered  photons.  Instead,  the  observed  photons  from  the  Co-57  source 
likely  are  those  that  scatter  multiple  times  inside  of  the  collimator.  By  bouncing 
back  and  forth  between  the  AIM-70  of  the  closed  chambers,  these  photons  maintain 
a  similar  trajectory  to  that  at  which  they  were  emitted,  but  with  lower  energies.  In 
this  way  the  ML-EM  algorithm  is  able  to  recreate  an  image  of  the  Co-57  source,  even 
though  the  energy  window  is  set  much  lower  than  122  keV. 

In  the  future,  this  effect  could  be  characterized,  leading  to  a  more  effective  system. 
Ideally,  the  full  energy  spectrum  of  detected  photons  would  be  used  to  reconstruct 
an  image  of  the  source.  Again,  with  a  larger  energy  window,  count  rates  increase. 
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Figure  44.  Images  of  the  Cd-109  button  source  with  background  using  Brown  correla¬ 
tion,  ML-EM  with  150  iterations,  and  ML-EM  with  500  iterations.  1,  5,  and  10  mask 
sequences  are  shown. 
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Returning  to  the  main  purpose  of  the  test,  the  use  of  multiple  mask  patterns 
successfully  reduced  noise  in  the  reconstructed  images.  The  inherent  noise  in  the 
correlation  images  is  greatly  reduced  with  the  use  of  multiple  masks.  With  the  ML- 
EM  method,  the  Cd-109  source  is  apparent  above  the  blur  from  the  Co-57  source. 
When  1  mask  was  used,  significant  noise  is  visible  in  regions  of  the  image  where  no 
source  was  present.  This  noise  is  reduced  when  more  masks  are  used.  The  Co-57  blur 
remains,  regardless  of  the  number  of  masks  used.  With  more  masks  the  size  of  the 
Cd-109  portion  of  the  image  is  honed  down. 

Using  the  SNR  analysis  methods  presented  in  Section  5.3,  the  SNR  and  SNRRMs 
were  found.  The  results  for  the  correlation  method  are  shown  in  Figure  45.  The 
SNRrms  increased  by  a  a  factor  of  3.04  when  ten  masks  were  used.  Notice  that 
the  discontinuities  between  2  and  3  mask  sequences  and  5  and  6  masks  sequences 
are  present  again,  as  they  were  in  the  SNR  plots  for  the  Co-57  images.  A  potential 
explanation  for  these  discontinuities  is  given  in  Section  5.5. 
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Figure  45.  Plot  of  the  SNR  versus  the  number  of  masks  used.  Correlation  method. 
Cd-109  button  source  with  background.  1  hour  total  measurement  time.  Error  bars 
are  shown. 
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The  ML-EM  SNR  plots  both  show  significant  increases  in  SNR.  SNR  increased 
by  a  factor  of  1.8  when  ten  masks  were  used.  The  discontinuity  between  the  5  and  6 
mask  SNR  results  is  possibly  due  to  the  explanation  given  in  Section  5.5.  ffowever, 
the  reason  for  the  SNR  jump  between  1  and  2  masks  is  unclear. 
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Figure  46.  Left:  SNR  versus  number  of  masks  used,  Cd-109  with  background,  1  hr  de¬ 
tection  time.  Right:  SNRrms  versus  number  of  masks  used,  Cd-109  with  background, 
1  hr  detection  time. 


This  source  configuration  is  somewhat  strange,  in  that  much  of  what  was  consid¬ 
ered  noise  originated  from  the  Co-57  source  located  within  the  field-of-view  of  the 
collimator/detector.  It  was  not  expected  that  the  Co-57  would  still  be  visible  in  the 
reconstructed  images.  The  increase  in  image  quality  is  clear  to  the  eye,  with  the 
Cd-109  source  becoming  sharper  when  more  masks  were  used.  For  a  future  test,  the 
source  of  background  radiation  should  be  placed  behind  the  detector  to  provide  a 
source  of  uniform  radiation. 

5.4.4  Co-57  Intermittent  Rod 

With  the  intermittent  rod,  about  100,000  events  were  recorded  by  the  detector  in 
a  given,  total  datarun.  Two  lead  bricks  were  placed  in  front  of  the  rod,  which  limited 
the  width  of  the  source  plane  to  about  7  cm.  This  was  done  to  remove  counts  that 
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would  arise  from  portions  of  the  source  outside  of  the  collimator’s  held-of-view.  The 
rod  has  a  width  of  2  mm  and  has  a  dotted  structure  along  its  length.  1  cm  long  lead 
pieces  are  wrapped  around  the  rod,  spaced  1  cm  from  each  other  down  the  length  of 
the  rod.  Moving  down  the  length  of  the  rod,  7-rays  should  be  emitted  for  1  cm,  then 
blocked  for  1  cm,  then  emitted  for  1  cm,  and  so  forth.  With  the  lead  bricks  in  place, 
only  four  emitting  segments  were  visible  to  the  detector.  These  four  segments  should 
be  apparent  in  the  reconstructed  images. 
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Figure  47.  Images  of  the  intermittent  rod  using  Brown  correlation,  ML-EM  with  150 
iterations,  and  ML-EM  with  800  iterations.  1,  5,  and  10  mask  sequences  are  shown. 
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The  correlation  methods  poorly  reconstruct  the  images  of  the  rod.  The  inherent 
noise  does  decrease  when  more  mask  pattens  are  used,  but  the  signal  portion  is  still  not 
very  good.  The  correlation  techniques  suffer  drastically  from  the  vignetting  effects. 
Most  coded-apertures  are  nowhere  near  as  deep  as  the  programmable  collimator. 
Most  of  the  effort  for  this  project  was  given  to  the  ML-EM  method.  Compensation 
for  the  vignetting  effects  in  the  correlation  methods  could  be  included  in  a  future 
test,  if  desired. 


Figure  48.  Left:  Picture  of  the  Rod.  Top-Right:  The  ML-EM  image  reconstruction 
using  1  mask,  zoomed.  Bottom-Right:  The  ML-EM  image  reconstruction  using  10 
masks,  zoomed. 


The  ML-EM  reconstructed  images  are  much  better.  With  1  mask  some  of  the 
segmentation  in  the  rod  is  visible.  With  10  masks  the  emitting  segments  of  the  rod 
are  clearly  visible.  The  10-mask  image  is  a  significant  improvement.  Notice  also  that 
the  maximum  value  of  the  colorscale  decreases  when  more  masks  are  used.  The  1- 
mask  image  is  peaked  in  an  unphysical  way.  The  emitting  segments  emit  uniformly 
across  their  lengths,  as  shown  in  the  10-mask  image. 
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5.4.5  Co-57  Rope 


The  Co-57  rope  was  imaged,  but  only  with  up  to  a  maximum  of  five  masks.  No 
shift  correction  was  performed  on  the  rope  images.  As  with  the  intermittent  rod,  two 
lead  bricks  were  placed  in  front  of  the  rope,  limiting  the  width  of  the  source  plane  to 
about  7  cm.  About  31,000  events  were  recorded  for  each  total  datarun,  within  the 
energy  window.  The  rope  lacks  the  intermittent  structure  of  the  rod.  The  width  of 
the  rope  is  2  mm.  Again,  the  correlation  method  did  not  perform  very  well.  Inherent 
noise  decreased  when  more  mask  patterns  were  used,  but  the  source  is  only  visible 
at  the  center  of  the  reconstructed  image  because  of  vignetting  effects.  The  ML-EM 
method  reconstructed  the  image  fairly  well. 


1  Mask  3  Masks  5  Masks 


Figure  49.  Zoomed  view  of  the  rope  after  1,  3,  and  5  masks,  ML-EM,  150  iterations. 


The  1-mask  ML-EM  reconstructed  image  is  blotted  in  an  unphysical  way.  The 
source  strength  should  be  uniform  along  the  length  of  the  rope.  This  blotting  is  less 
noticeable  in  the  5-mask  image.  In  some  ways  the  3-mask  image  is  better  than  the 
5-mask  image,  especially  with  regard  to  the  lower  half  of  the  rope.  A  zoomed  version 
of  this  lower  portion  is  given  in  Figure  49.  Unfortunately,  this  source  was  only  imaged 
with  up  to  five  masks.  If  ten  masks  had  been  used,  the  final  image  would  have  been 
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much  better.  In  the  case  of  the  rod,  most  of  the  improvement  occurred  between  5 
and  10  mask  sequences. 


Figure  50.  Images  of  the  rope  using  Brown  correlation,  ML-EM  with  150  iterations, 
and  ML-EM  with  800  iterations.  1,  3,  and  5  mask  sequences  are  shown. 


The  rope  takes  up  more  of  the  source  plane  than  did  the  rod,  resulting  in  more 
noise.  The  coded-aperture  multiplex  advantage  decreases  as  a  larger  portion  of  the 
held  of  view  is  taken  up  by  the  source.  Nevertheless  the  image  did  improve  when 
more  masks  were  used. 


5.5  Discontinuities  in  the  SNR  from  Correlation  Decoding 


16 

14- 

12' 

CO 

of  10- 

z 

in 

8- 

6 

4- 


4  6  8  10 

Number  of  Masks 


Figure  51.  Left:  SNRrms  f°r  the  Co-57  button  source,  as  a  function  of  number  of  masks 
used,  correlation  decoding.  Right:  SNRrms  f°r  the  Cd-109  source  with  background, 
as  a  function  of  number  of  masks  used,  correlation  decoding. 


Regarding  images  reconstructed  with  correlation  methods,  significant  discontinu¬ 
ities  are  present  in  the  plots  of  SNRrms  versus  the  number  of  masks  used.  For  both 
the  Co-57  source  and  the  Cd-109  source,  these  discontinuities  occur  between  2  and  3 
masks  and  between  5  and  6  masks.  First  it  was  thought  that  this  effect  was  caused 
by  the  way  in  which  the  signal  circle  was  chosen.  The  circle  which  cuts  out  the  signal 
portion  of  the  image  is  chosen  qualitatively.  If  this  circle  is  too  small,  then  the  varia¬ 
tions  in  the  signal  shape  might  cause  a  piece  of  the  signal  to  lie  outside  of  the  circle, 
meaning  that  some  of  the  signal  would  be  counted  as  noise.  To  test  for  this  effect,  the 
radius  of  the  circle  was  increased.  The  SNRrms  results  were  about  the  same,  though 
shifted  downward  due  to  more  noise  being  included  in  the  signal  list.  Ultimately,  the 
actual  sizes  of  the  sources  were  used  for  the  signal  circle. 

Another  possibility  considered  was  that  masks  3  and  6,  chosen  at  random,  just 
so  happened  to  have  good  correlation  properties.  Recall  from  Section  2.1  that  the 
term  A-kG  should  approach  a  h-function  as  nearly  as  is  possible,  where  A  is  the  mask 
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Figure  52.  The  average  noise  in  A*G  versus  the  number  of  masks  used.  Here  lower 
values  mean  that  A  *  G  is  closer  to  a  (5-function 

pattern,  G  is  the  Brown  matched  filter,  and  *  is  the  2D  cross-correlation  operator. 
The  correlation  of  each  mask  pattern  with  its  matched  filter  was  calculated.  When 
a  Brown  filter  is  used,  the  central  spike  of  A  *  G  always  will  be  zero.  As  the  noise 
in  the  rest  of  the  matrix  decreases,  becoming  more  negative,  the  term  more  closely 
approximates  a  5-function.  These  results  are  shown  in  Figure  52.  Certainly  some 
of  the  mask  patterns  have  better  correlation  properties  than  others.  This  can  be 
considered  in  future  work.  However,  the  shape  of  the  curve  in  Figure  52  does  not 
explain  the  discontinuities  in  Figure  51. 


0.515 


0.51  -  • 

CD 

-f  0.505 

CD 

“0.5  - 

©  \ 

Q. 

g  0.495 

Q. 

O 

©  0.49  •  | 

o>  0.485 
> 

< 

0.48  i 

0.475  ■ - ■ - * - - - 

0  2  4  6  8  10 

Number  of  Masks  in  Sequence 


Figure  53.  Left:  The  average  correlation  of  a  given  mask  and  all  other  masks  versus 
the  number  of  masks  used.  Right:  The  fraction  of  time  in  which  the  average  chamber 
is  open  to  y-rays  versus  the  number  of  masks  used. 


In  fact,  the  data  presented  in  Figure  53  are  needed  to  explain  the  discontinuities. 
The  left-hand  plot  shows  the  average  cross-correlation  between  one  mask  and  all  of 
the  other  masks  in  the  sequence.  The  right-hand  plot  shows  the  fraction  of  time  in 
which  the  average  chamber  is  open,  as  a  function  of  the  number  of  masks  used.  Notice 
that  for  sequences  greater  than  1  mask,  the  two  plots  are  identical  but  for  ordinate 
scaling.  This  is  not  surprising,  given  the  relationship  that  was  established  between 
mask  independence  and  the  fractional  time  that  each  chamber  is  open,  as  explained 
in  Section  3.3.  It  was  shown  that  the  least  correlated  sequence,  or  the  sequence  in 
which  the  mask  patterns  were  most  independent,  is  that  in  which  each  chamber  is 
“on”  and  “off”  for  equal  durations.  To  have  non-trivial  mask  patterns,  the  “on”  and 
“off”  durations  cannot  be  guaranteed  as  equal  except  at  the  end  sequence  with  the 
maximum  number  of  masks. 


Figure  54.  The  fraction  of  time  in  which  the  average  chamber  is  open  was  subtracted 
from  the  ideal  fractional  duration,  0.5.  The  negative  of  the  absolute  value  was  plotted. 


Perhaps  the  most  enlightening  version  of  the  data  is  given  by  Figure  54.  In 
this  case,  the  fractional  time  in  which  the  average  collimator  chamber  was  “on”  was 
subtracted  from  the  ideal  fractional  duration  of  0.5.  The  negative  absolute  value  of 
this  difference  was  then  taken.  Notice  how  similar  Figure  54  is  to  Figure  51.  Clearly 


the  fractional  duration  for  which  the  average  chamber  is  “on”  implicitly  drives  the 
degree  to  which  SNR  increases. 

5.6  MCST  Uncollimated  Forward  Problem 

5.6.1  Results 

The  different  metal  samples  were  examined  for  different  periods  of  time.  It  was 
found  empirically  that  about  24  hours  of  data  were  required  to  get  significant  Gaus- 
sians  in  each  of  the  subpixels,  given  the  setup  used.  In  the  first  test  to  be  evaluated, 
a  copper  strip  was  secured  to  the  fishing  line  and  data  were  taken  for  24  hours. 

An  energy  spectrum  for  each  of  the  detector  pixels  was  required  from  the  detector 
software.  In  one  mode  the  software  will  export  all  of  the  events  that  it  detects, 
between  a  certain  energy  window,  into  a  text  hie.  Originally  the  window  was  set  so 
that  events  with  energies  between  60  keV  and  80  keV  were  recorded.  Later  this  was 
changed  to  accept  events  between  53  keV  and  83  keV,  so  as  to  allow  the  Gaussian 
fit  routine  to  have  a  better  baseline.  Each  line  of  the  hie  corresponds  to  a  given 
detector  event  and  gives  energy,  event  pixel,  etc.  A  parser  was  written  in  Matlab® 
to  transform  the  data  from  the  raw  events  into  energy  histograms  for  each  detector 
pixel.  These  data  were  stored  as  a  three-dimensional  mat-hle. 

Matlab®  is  equipped  with  a  built-in  function  for  fitting  a  Gaussian  distribution 
to  data  points  using  a  least-squared  algorithm.  The  function  does  not  allow  for  an 
offset  in  the  data,  however.  Another  script  was  written  which  hrst  imports  the  mat- 
hle  containing  the  energy  spectra.  For  each  spectrum,  the  script  finds  the  average  of 
the  hrst  ten  data  points  as  well  as  the  average  of  the  last  ten  data  points.  A  line  is 
made  between  the  two  averages,  and  this  line  is  subtracted  from  all  of  the  data.  Next, 
the  Gaussian  ht  function  is  performed  for  all  spectra.  The  line  is  then  added  back 
to  both  the  ht  and  the  raw  data.  The  Gaussian  function  outputs  the  Full  Width  at 
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The  energy  spectrum  at  detector  pixel  (8,8).  100-bin  histogram 


The  peak  energy  of  the  gaussianfitat  each  detector  pixel 
11101 9_MC ST_C d  1 09_C u  -  colorbar  in  keV 


2  4  6  8  10  12  14  16 

Detector  Column 


The  counts  under  the  Gaussian  fit,  colorbar  in  counts 
1 1 1019_MCST_Cd109_Cu 


Figure  55.  Selected  data  from  a  test  using  a  single  copper  sample.  Compare  with 
the  model  predictions  shown  in  Figure  16.  Top:  The  energy  spectrum  from  a  single 
detector  pixel,  with  its  Gaussian  fit.  Left:  The  energy  of  the  peak  of  the  fitted  Gaussian 
in  each  detector  pixel.  Right:  The  number  of  counts  under  each  fitted  Gaussian  for 
each  detector  pixel. 
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Half  Maximum  (FWHM)  and  the  Root  Mean  Squared  (RMS)  error  of  the  fit,  among 
other  things.  The  number  of  counts  for  each  pixel  is  an  important  number  for  MOST. 
The  Gaussian  fit  is  integrated  using  an  adaptive  quadrature  method,  also  built-in  to 
MATLAB®  .  This  integration  gives  the  number  of  counts  attributable  to  Compton 
events  in  the  sample. 

The  pixels  that  lie  on  the  corners  of  the  detector  will  receive  significantly  lower 
counts  than  the  rest  of  the  pixels,  because  they  are  physically  smaller.  For  the  exper¬ 
iments  conducted,  these  corner  pixels  typically  received  too  few  counts  for  a  Gaussian 
fit  to  be  reliable.  Thus,  the  data  and  the  fits  from  these  pixels  were  discarded.  With 
higher  count  rates  these  pixels  could  be  utilized.  However,  the  model  must  first  take 
their  lower  effective  efficiencies  into  account. 

Raw  and  Gauss  lit,  pixel  (4,1).  An  Example  "Comer  Pixel",  whose  fit  should  be  ignored 


190CT1 1  ,MCST,Cd109,CuStrip  -  24  hour  run 


Figure  56.  One  of  the  corner  pixels,  for  the  same  data  shown  in  Figure  55.  Notice  the 
low  number  of  counts. 

In  the  results  shown  in  Figure  55,  where  a  copper  sample  was  examined,  the 
average  FWHM  of  the  Gaussian  fits  was  found  to  be  4.91  keV.  The  average  RMS 
error  of  the  fits  was  5.36  counts.  The  average  number  of  counts  into  each  pixel  was 
366.  Something  to  consider  is  that  the  samples  used  are  not  as  small  as  the  potential 


92 


resolution  of  the  system.  FHWM  values  add  in  quadrature.  That  is, 


Rtot  =  sjRl  +  Rl  (32) 

where  R  stands  for  resolution,  given  here  by  the  FWHM.  The  detector  has  an  energy 
resolution  of  about  1.4  keV  at  88  keV,  as  was  found  in  a  test.  The  size  of  the  metal 
could  be  reduced  substantially  before  the  limits  of  the  detector  are  reached. 

In  another  test  the  copper  sample  was  cut  down  to  a  quarter  of  its  former  area. 
All  other  parameters  were  held  fixed.  It  was  thought  that  the  FWHM  would  decrease. 
However,  the  average  FWHM  was  calculated  as  5.74  keV,  with  a  fit  RMS  error  of 
5.44  counts.  The  average  number  of  counts  was  288.  And  so,  by  quartering  the  size  of 
the  copper  strip,  the  resolution  did  not  change  much,  although  the  number  of  counts 
received  by  the  detector  did  change.  This  means  that  the  resolution  of  an  MOST 
system  could  be  very  small,  if  perfected. 

In  a  further  experiment  aluminum  was  used  instead  of  copper.  Aluminum  has 
a  significantly  lower  Z  than  copper  (13  versus  29).  One  would  expect  the  counts 
received  from  the  aluminum  to  be  significantly  lower  than  from  copper.  Two  tests 
were  conducted,  both  over  24  hours,  with  the  same  geometry.  The  resolution  and 
RMS  error  between  the  two  tests  was  comparable.  The  average  number  of  counts 
received  from  the  aluminum  sample  was  found  to  be  272  which  is  substantially  lower 
than  the  366  counts  received  by  the  average  detector  pixel  in  the  copper  test.  These 
results  indicate  that  the  system  could  successfully  differentiate  between  materials 
with  different  electron  densities. 

A  CuCu  and  an  AlCuCu  (Figure  20)  phantom  were  tested  using  the  system,  as 
described  in  section  4.3.  The  results  from  each  are  comparable,  with  those  from  the 
AlCuCu  phantom  shown  in  Figure  57.  The  built-in  Matlab®  routine  has  the  ability 
to  fit  multiple  Gaussian  distributions  to  data.  However,  it  was  wholly  unsuccessful  at 
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fitting  multiple  Gaussian  distributions  to  these  data.  The  Gaussian  distributions  are 
simply  too  close  together,  in  energy,  for  that  routine  to  work.  More  robust  methods 
are  necessary.  The  PWLS  algorithm  is  one  choice,  although  ML-EM  may  be  usable 
as  well. 


Detector  Pixel  (6,9),  25  hour  run,  Cd-109 
1 11021  data,  AlCuCu  Phantom 
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Figure  57.  Spectra  received  from  the  AlCuCu  phantom,  taken  from  two  different 
detector  pixels. 


5.6.2  Discussion 

The  results  found  from  the  basic  MCST  setup  are  in  good  agreement  with  the 
MCST  forward  model  predictions.  A  change  in  detected  photon  energy  across  the  face 
of  the  detector  crystal  is  certainly  apparent.  The  coded- aperture  results  give  some 
insight  into  how  the  programmable  collimator  might  help  with  MCST.  In  MCST, 
photons  are  scattered  within  the  sample  and  the  scattered  photons  are  then  detected. 
Effectively  the  sample  can  be  thought  of  as  an  extended  source,  similar  to  the  rod  and 
the  rope  as  presented  in  Section  5.4.5.  Just  as  the  programmable  collimator  improved 
the  images  of  the  rod  and  rope,  so  it  would  improve  the  images  of  an  MCST  source. 

In  the  coded-aperture  experiments,  photon  energy  is  removed  from  consideration. 
The  energy  of  the  photons  emitted  from  the  radionuclide  was  known,  and  the  detector 
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was  set  only  to  record  photons  of  that  energy  (plus  or  minus  a  few  keV).  In  this  way, 
Compton  scattered  photons  were  removed  from  consideration,  resulting  in  less  noisy 
images.  In  MOST,  however,  the  photon  energy  is  needed  in  the  reconstruction  process. 
This  requirement  introduces  another  dimension  of  complexity  that  was  not  considered 
in  the  coded-aperture  experiments. 

The  HPGe  strip  detector  used  in  this  study  is  uniquely  suited  to  handle  MOST. 
HPGe  detectors  have  very  good  energy  resolution  [37].  The  pixel  size  and  number  of 
pixels  is  comparable  to  detectors  made  of  other  materials.  Both  position  sensitivity 
and  precise  energy  sensitivity  are  required  for  MOST.  The  one  disadvantage  of  the 
HPGe  strip  detector  is  its  low  efficiency.  Perhaps  the  efficiency  could  be  increased 
using  Compton  rescue  [67]. 

Evans  discovered  that  for  three-dimensional  MCST,  some  sort  of  physical  con¬ 
straint  was  required,  typically  taking  the  form  of  collimation  [25].  Without  collima- 
tion,  the  ill-conditioned  matrices  would  not  converge  during  the  image  reconstruction 
process.  With  collimation,  the  parameters  could  be  sufficiently  bounded  to  allow  for 
convergence.  This  is  where  the  programmable  collimator  applies  to  MCST.  The  term 
multiplexing  has  been  used  throughout  this  document.  When  applied  to  MCST,  mul¬ 
tiplexing  simply  means  that  each  detector  pixel  receives  information  about  more  than 
one  pixel  in  the  sample.  The  sample  is  discretized,  and  if  any  of  the  detector  pixels 
detect  scattered  photons  from  more  than  one  of  the  sample  pixels,  then  multiplexing 
is  occurring. 

A  balance  must  be  struck  between  multiplexing  and  the  reconstruction  program’s 
ability  to  converge  on  a  solution.  More  multiplexing  means  lower  detection  time.  In 
medical  imaging,  more  multiplexing  means  a  lower  patient  dose  rate  of  ionizing  radia¬ 
tion.  In  materials  investigation,  more  multiplexing  could  mean  real-time  images  with 
good  quality.  On  the  other  hand,  as  mentioned  above,  if  the  multiplexing  becomes 
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Figure  58.  The  level  of  multiplexing  is  determined  by  the  sample-collimator  and 
collimator-detector  distances.  Top:  little  to  no  multiplexing.  Middle:  medium  multi¬ 
plexing.  Bottom:  High  multiplexing. 
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too  great,  the  ill-conditioned  matrices  will  not  converge  during  the  image  reconstruc¬ 
tion  process.  Or,  if  convergence  does  occur,  the  images  might  be  noisy  and  unusable 
if  there  is  too  much  multiplexing. 

The  programmable  collimator  could  be  used  to  limit  the  problem  while  still  al¬ 
lowing  for  multiplexing.  In  the  coded-aperture  experiments,  the  collimator  depth 
was  a  detriment.  In  MOST,  the  collimator’s  depth  could  be  used  as  an  advantage, 
providing  geometric  constraints.  By  adjusting  the  sample-collimator  and  collimator- 
detector  distances,  the  level  of  multiplexing  is  adjusted.  The  multiplexing  can  be 
tinkered  with  by  the  experimenter  until  the  optimal  setup  is  found.  Then,  by  using 
the  programmable  collimator,  the  overall  SNR  should  increase  because  of  the  mul¬ 
tiple,  semi-independent  views  afforded  to  the  detector.  Just  as  the  images  of  the 
extended  sources  in  the  coded-aperture  experiments  improved,  so  would  the  MOST 
images  improve  when  multiple  mask  patterns  are  used.  The  multiple  masks  would 
further  chop  up  the  isogonic  arcs  found  from  the  Compton  formula  (Figure  7). 

5.7  Programmable  Collimator  Physical  Performance 

The  use  of  the  programmable  collimator  was  hampered  by  several  factors  through¬ 
out  its  life.  First,  the  square  chambers  were  a  bit  more  difficult  to  machine  than  cir¬ 
cular  chambers  would  have  been.  It  was  impossible  for  the  wire  EDM  to  cut  chambers 
with  perfectly  square  cross  sections  into  the  grid  block.  The  corners  were  rounded. 
The  square  PEEK  plugs  had  to  be  hied  and  broached  individually.  In  effect,  each 
plug  was  custom-made  for  its  particular  chamber.  Although  the  square  chambers  are 
better  for  calculations  and  allow  more  throughput,  circular  holes  might  be  considered 
in  future  designs  to  decrease  construction  time. 

Secondly,  the  plugs  had  to  fit  tightly  into  the  chambers  to  prevent  AIM-70  from 
leaking  out  of  the  front  of  the  device.  Often  they  were  so  tight  that  the  plugs  would 
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break  when  being  pulled.  There  was  a  danger  of  pulling  a  given  plug  all  the  way 
out  of  its  chamber,  as  described  in  Section  4.4.  To  prevent  spillage  of  AIM-70,  the 
collimator  was  evacuated  of  AIM-70  before  the  mask  pattern  was  reconfigured.  The 
collimator  was  then  refilled  with  AIM-70.  The  entire  process  of  changing  the  mask 
configuration  took  from  1-2  hours.  Such  a  lengthy  procedure  is  tolerable  for  a  proof- 
of-concept  experiment,  but  a  better  system  is  desirable. 

Finally,  portions  of  the  grid  block  cracked  and  chipped  away  during  the  pro¬ 
grammable  collimator’s  use.  This  was  probably  clue  to  the  stress  placed  on  the  grid 
walls  when  the  AIM-70  cooled  and  expanded.  Potential  solutions  are  given  in  Section 


VI.  Conclusion  and  Path  Forward 


6.1  Summary 

The  first  fully  reconfigurable  programmable  collimator  using  a  liquid  attenuator 
was  designed,  built,  and  tested.  The  programmable  collimator  was  considered  in  two 
modes,  as  a  coded-aperture  and  as  a  problem  constraining  device  for  MOST.  For  the 
coded-aperture  application,  an  MCNP  simulation  was  conducted  in  order  to  better 
understand  the  problem.  A  M  ATLAB®  script  was  written  to  find  the  sequence  of  least- 
correlated  random  mask  patterns,  so  as  to  have  the  most  independent  views  possible 
of  the  object  plane.  A  ray-tracing  model  was  created  in  Matlab®  to  characterize 
the  forward  problem,  providing  the  necessary  information  to  the  ML-EM  algorithm. 
Vignetting  effects,  the  shape  of  the  detector  crystal  and  dead  subpixels  were  factored 
in  to  the  forward  problem,  increasing  accuracy.  An  experiment  was  constructed  to 
test  the  device  as  a  coded-aperture. 

Three  decorrelation  techniques  were  evaluated  using  the  experimental  data.  The 
decorrelation  techniques  were  useful  because  of  their  simplicity  and  computational 
speed.  The  ML-EM  method  was  also  evaluated.  It  generally  performed  better  than 
the  decorrelation  techniques  with  respect  to  SNR,  though  the  greater  number  of 
tweak-able  parameters  makes  the  method  significantly  more  complex.  ML-EM  is 
also  much  more  computationally  expensive  than  the  correlation  methods.  Regardless 
of  the  method,  the  programmable  collimator  is  demonstrably  effective  as  a  coded- 
aperture.  A  sequence  of  ten  mask  patterns  was  used,  and  the  SNR  increased  with  the 
number  of  mask  patterns  used.  This  promising  result  suggests  that  the  programmable 
collimator  could  significantly  improve  coded-aperture  technology.  A  sequence  with 
twenty  or  more  mask  patterns  should  be  tested  in  the  future  in  order  to  determine  at 
how  many  masks  the  SNR  will  peak  over  a  constrained,  total  detection  time. 
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A  model  was  developed  for  MOST,  based  on  the  Compton  equation  and  the  Klein- 
Nishina  formula.  An  experiment  was  conducted  in  which  small  samples  were  placed 
in  front  of  a  mono-energetic  source  and  the  scattered  photons  were  measured.  The 
experimental  results  suggest  that  Compton  scatter  tomography  is  possible  with  the 
PHDs  DSSD.  A  robust  method  like  PWLS  would  have  to  be  used  in  order  to  unfold 
the  individual  responses  from  all  of  the  sample  voxels.  In  the  past,  MCST  has  failed 
because  it  is  too  highly  multiplexed  and  thus  becomes  ill-conditioned  and  unsolvable. 
The  programmable  collimator  may  be  able  to  physically  constrain  the  problem,  al¬ 
lowing  convergence  on  a  solution.  The  problem  would  still  be  multiplexed  to  some 
extent,  but  not  so  multiplexed  as  to  be  unsolvable.  By  having  some  multiplexing, 
performance  would  increase  relative  to  the  highly  collimated  systems  used  today. 

6.2  Practical  Uses 

Please  see  Sections  2.5  and  2.6.5  for  background  information. 

6.2.1  MCST 

The  programmable  collimator  should  be  a  viable  constraining  device  for  MCST. 
MCST  has  several  potential  uses.  In  medical  imaging,  Compton  scatter  dominates  as 
the  mode  of  interaction  in  the  patient’s  body  [7].  Most  systems  ignore  and  attempt 
to  filter  out  the  photons  that  were  Compton  scattered  within  the  subject.  If  those 
events  were  utilized,  image  fidelity  would  increase  and/or  dose  to  the  patient  would 
decrease.  3-dimensional  images  could  be  made  more  readily  with  the  additional  in¬ 
formation  gained  from  Compton  scatter  events.  Likely,  a  conventional  system  would 
be  augmented  with  a  programmable  collimator  detector. 

Another  potential  system  would  measure  Compton  scattered  photons  exclusively. 
Such  a  system  would  be  useful  for  industrial  cases  in  which  access  to  only  one  side 
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of  a  sample  is  available.  In  this  case  the  physical  collimation  from  the  programmable 
collimator  would  be  a  necessity.  Without  collimation,  the  problem  is  too  multiplexed 
to  solve.  The  programmable  collimator  would  help  to  decrease  detection  time  or 
source  intensity. 

6.2.2  Unmanned  robots  or  UAVs 

A  detection  system  utilizing  programmable  collimation  would  be  useful  for  UAV 
and  unmanned  robot  applications.  The  unmanned  vehicles  could  be  used  to  survey 
large  areas  after  a  nuclear  explosion  or  subsequent  to  wide-area  dispersion  of  nuclear 
material  following  a  dirty  bomb  attack  or  a  reactor  accident.  Inexpensive  position 
sensitivity  would  help  to  limit  the  time  required  to  map  the  contamination  over  a  large 
area.  Lower  detection  times  are  especially  enticing  given  battery  life  constraints  in 
remotely  controlled  devices.  Unmanned  vehicles  could  also  be  used  to  sniff  out  nuclear 
weapons  before  they  are  used,  if  other  intelligence  provides  a  general  location  for  the 
nuclear  device.  Inexpensive  position  sensitivity  would  help  to  find  the  nuclear  material 
more  quickly.  Also,  collection  of  radioactive  effluents  is  an  important  component  of 
several  treaty-monitoring  organizations.  A  position-sensitive  system  would  help  to 
guide  members  of  these  organizations  toward  the  proper  locations  where  effluent 
content  is  highest. 

6.2.3  Space  Surveillance 

Programmable  collimation  would  not  be  advantageous  for  space-based  detection 
of  nuclear  events  on  the  ground.  The  duration  of  7-ray  output  from  a  nuclear  weapon 
is  too  short  for  a  changing  collimator  to  be  of  any  real  use.  However,  the  idea 
of  using  multiple  masks  is  still  valid.  Several  small  satellites  could  be  equipped 
with  collimators  that  have  different  mask  patterns.  The  mask  patterns  would  be 
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chosen  such  that  each  satellite  would  receive  the  most  unique  information  possible. 
Knowing  the  location  of  each  satellite  precisely,  the  images  from  each  satellite  could 
be  combined  to  form  a  better  image  than  could  be  found  from  a  single  satellite. 

6.2.4  Space-based  Gamma-Ray  Astronomy 

Programmable  collimation  most  certainly  could  be  useful  for  7-ray  astronomy. 
Gamma  Ray  Bursts  (GRBs)  have  durations  on  the  order  of  minutes,  meaning  that 
many  mask  patterns  could  be  used  as  the  event  occurs.  SNR  and  image  resolution 
of  the  GRBs  would  thus  increase  significantly  relative  to  current  imaging  systems. 
For  many  years  coded-apertures  have  been  used  aboard  satellites  for  astronomy,  suc¬ 
cessfully.  Programmable  coded-apertures  are  a  logical  next-step.  An  array  of  smaller 
satellites  equipped  with  programmable  coded-apertures  could  outperform  the  mas¬ 
sive,  bulky  satellites  used  currently.  The  reconfigurability  of  the  collimator  would 
eliminate  the  need  for  the  satellite  rotation. 

6.2.5  Medical  Nuclear  Imaging 

Coded-apertures  have  been  effectively  demonstrated  in  Single  Photon  Emission 
Computed  Tomography  (SPECT)  and  Positron  Emission  Tomography  (PET).  Pro¬ 
grammable  coded-apertures  represent  a  logical  progression  in  that  branch  of  imaging. 
Regardless  of  the  coded-aperture  setup  used,  programmable  collimation  can  increase 
the  image  resolution.  In  medical  situations  there  is  adequate  time  for  multiple  mask 
patterns  to  be  used.  The  ultimate  result  of  programmable  collimation  would  be  better 
images  or  lower  dosage  given  to  patients. 
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6.2.6  Position  Sensitivity  from  a  Single  Detector  Crystal 

Perhaps  the  most  enticing  application  for  the  programmable  collimator  is  that 
it  could  give  position  sensitivity  to  a  single  detector.  If  this  could  be  achieved,  the 
price  for  a  position-sensitive  detector  could  be  reduced  substantially.  An  ideal  setup 
would  be  to  combine  a  Nal  detector  and  a  programmable  collimator.  The  Nal  detec¬ 
tor  has  relatively  good  energy  resolution.  For  many  applications,  cost  is  the  limiting 
factor  with  regard  to  position  sensitive  detectors.  An  inexpensive,  position  and  en¬ 
ergy  sensitive  detector  would  be  helpful  in  all  of  the  applications  listed  above.  For 
the  space-based  astronomy  application,  Nal  detectors  also  have  much  lower  power 
requirements  than  HPGe  strip  detectors  and  are  more  resistant  to  vibration. 

The  ease  with  which  the  programmable  collimator  was  modeled  suggests  that 
position  sensitivity  could  be  given  to  single  detector.  Only  iterative  methods  could 
be  used  for  that  type  of  problem. 

6.3  Future  Work 

Throughout  this  project,  several  areas  have  been  highlighted  which  would  benefit 
from  further  investigation.  Based  on  the  current  setup,  one  could: 

•  Characterize  the  efficiency  of  each  subpixel  in  the  HPGe  strip  detector.  In 
this  work,  all  subpixels  were  assumed  to  have  equal  efficiencies.  This  is  not 
the  case  clue  to  constructional  variations,  varying  pixels  sizes,  and  nonuniform 
charge  movements.  A  characterization  experiment  would  require  that  a  highly 
collimated  source  be  precisely  moved  such  that  one  subpixel  is  irradiated  at  a 
time  [75].  Comparison  of  counts  would  provide  efficiency. 

•  Incorporate  more  effects  into  the  forward  model.  In  this  work  Compton  scat¬ 
tering  within  the  programmable  collimator  was  ignored,  as  well  as  attenuation 


103 


through  the  PEEK  plugs.  If  the  ML-EM  algorithm  is  supplied  with  a  more  ac¬ 
curate  model  that  includes  those  factors,  and  others,  it  would  provide  a  better 
estimation  of  the  source  plane. 

•  Use  more  active  sources  or  longer  detection  times  so  as  to  have  better  counting 
statistics.  In  this  research,  data  and  collection  times  were  limited.  An  increase 
in  either  of  those  areas  would  help  to  rule  out  poor  counting  statistics  as  a  cause 
for  artifacts  in  the  reconstructed  images.  Counting  statistics  were  especially 
poor  in  the  MOST  testing. 

•  Perform  an  experiment  to  determine  the  limit  on  the  number  of  mask  patterns 
used,  given  a  constant  total  detection  time.  In  this  work,  the  image  quality 
increased  as  more  mask  patterns  were  used.  However,  if  the  overall  detection 
time  is  limited,  then  each  mask  is  used  for  a  shorter  duration  as  the  number  of 
masks  increases.  Eventually,  if  too  many  masks  were  used,  there  would  be  too 
little  time  to  formulate  an  image  through  each  mask.  Thus,  the  performance 
increase  has  a  peak  with  respect  to  the  number  of  mask  patterns  used. 

•  Image  7-rays  of  higher  energies.  A  dynamic  collimator  was  considered  necessary 
for  high  energy  7-rays  that  are  able  to  penetrate  the  “off”  portions  of  the 
collimator  [16].  Multiple  mask  patterns  should  effectively  overcome  the  noise 
introduces  because  of  penetration  through  the  AIM-70,  but  that  experiment 
needs  to  be  performed. 

•  Determine  how  to  choose  mask  patterns  on  the  fly.  A  computer  routine  could 
be  developed  to  use  incoming  data  in  order  to  make  an  educated  decision  on 
the  remaining  mask  patterns.  The  collimator  could  start  with  a  high  pinhole 
density,  allowing  maximum  throughput.  Then,  as  a  rough  source  position  is 
determined,  pinhole  density  would  decrease  so  as  to  increase  resolution.  In  this 
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way  a  sharp  image  could  be  found  in  minimal  time  even  when  viewing  a  wide 
area. 

Changing  the  experimental  configuration,  one  could: 

•  Place  the  programmable  collimator  in  front  of  a  single,  Nal  detector.  It  may 
be  best  to  perform  simulations  first,  in  order  to  understand  the  best  setup  that 
would  provide  position  sensitivity  to  the  single  detector. 

•  Use  a  setup  similar  to  that  presented  in  Section  4.3,  testing  the  programmable 
collimator  as  a  physical  constraint  on  the  MOST  problem.  A  forward  model 
would  have  to  be  developed  first,  analogous  to  that  presented  in  Section  3.1. 

•  Image  a  three-dimensional  source.  In  the  forward  model  for  ML-EM,  the  source 
region  would  be  divided  into  3D  voxels,  rather  than  2D  pixels.  Multiple  mask 
patterns  and  changes  in  magnification  would  allow  a  3D  source  to  be  imaged 
without  moving  the  collimator-detector  assembly.  In  other  words  tomography 
could  be  performed  without  rotating  the  detector  around  the  sample. 

6.4  The  Future  of  Programmable  Collimation 

The  next  design  would  be  similar  to  the  current  design,  with  a  few  notable  dif¬ 
ferences.  The  collimator  would  have  a  greater  number  of  smaller  holes,  with  each 
chamber  being  about  2  mm  wide  so  as  to  better  conform  to  current  coded-aperture 
designs.  With  more  chambers,  the  URA  and  MURA  mask  patterns  could  be  at¬ 
tempted,  although  it  is  unclear  whether  or  not  those  classes  of  mask  patterns  would 
be  optimal  for  ML-EM.  It  may  be  that  a  better  class  of  mask  patterns  exists  for  ML- 
EM.  Secondly,  the  collimator  could  be  made  with  less  depth,  reducing  vignetting. 
Thirdly,  the  holes  could  be  made  round,  to  substantially  reduce  construction  time. 
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High  temperature  O-rings  could  be  placed  around  the  plugs  to  prevent  leakage  and 
to  allow  easy  movement  of  the  plugs.  Fourthly,  screw  patterns  could  be  tapped  into 
each  plug,  so  that  the  actuator  rod  could  screw  into  the  plugs  to  move  them.  This 
might  prevent  the  plugs  from  breaking,  though  stripping  might  be  an  issue. 

A  different  liquid  attenuator  should  be  considered  prior  to  the  construction  of  the 
next  device.  Mercury  is  ideal.  It  has  a  high  Z  (80),  making  it  a  good  7  attenuator. 
Plus,  mercury  is  already  a  liquid  at  room  temperature,  removing  the  need  for  heat¬ 
ing  apparatuses.  Unfortunately,  mercury  is  difficult  to  procure,  due  to  its  toxicity. 
Other  eutectics  or  near-eutectics  could  be  a  possibility.  A  few  are  shown  in  Table  7. 
The  lack  of  toxicity  would  be  helpful,  in  order  to  eliminate  the  need  for  a  hood  or 
breathing  equipment.  A  lower  melting  point  would  be  useful.  A  material  that  does 
not  expand  upon  solidification  is  desired,  in  order  to  prevent  long-term  damage  to 
the  programmable  collimator. 


Table  7.  Potential  Candidates  for  the  liquid  7  attenuator  [13,  6,  76], 


Name 

Melting  Point 

[C] 

Expansion  After  Solidifying 
[inch  per  inch] 

Toxic? 

AIM-70 

70 

+0.005 

Fumes, Ingest 

Mercury 

-38.83 

- 

Yes 

AIM-47 

47 

+0.0002 

Fumes, Ingest 

AIM-58 

58 

+0.0002 

Ingest 

The  ultimate  goal  for  the  programmable  aperture  is  a  device  with  many  more 
smaller  chambers  that  is  fully  computer-automated.  Toward  that  goal  an  experimen¬ 
tal  system  could  be  built  where  pipettes  full  of  the  liquid  attenuator  are  used.  The 
liquid  attenuator  could  manually  be  squeezed  in  and  out  of  the  chambers.  Holes 
perpendicular  to  the  chambers  would  be  drilled  in  the  PEEK  (or  another  lightly  at- 
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tenuating,  machinable  plastic)  on  each  side  of  the  chamber  block,  forming  two  mani¬ 
folds.  The  manifolds  would  allow  the  liquid  to  be  pumped  either  into  the  chambers 
or  out  and  away  from  the  grid  block.  Eventually  a  pump  system  would  be  needed 
to  fully  automate  the  system.  The  attenuator  could  be  a  molten  metal,  as  in  this 
study,  a  liquid  metal  like  mercury  or  something  else  entirely.  The  chambers  would 
be  fabricated  to  a  small  enough  diameter  so  that  the  capillary  action  in  the  liquid 
attenuator  would  match  gravity,  ensuring  that  each  chamber  is  filled  vertically.  This 
would  allow  for  any  possible  orientation  of  the  device.  Either  positive  air  pressure  or 
a  non-attenuating  liquid  would  be  needed  on  both  sides  to  hold  the  liquid  attenuator 
in  the  chamber.  Low  cohesion  between  the  liquid  attenuator  and  the  other  parts  is 
desired  so  that  when  a  given  chamber  is  evacuated  no  liquid  attenuator  remains. 
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Attenuating 
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1 

■ 
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Figure  59.  A  future  programmable  collimator.  A  pumping  system  would  fill  or  evacuate 
each  chamber  with  an  attenuating  liquid.  Plastic  manifolds  on  each  end  of  the  grid 
block  would  guide  the  liquid  from  the  pumping  system  into  the  chambers. 
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A  device  that  can  quickly  switch  to  new  mask  configurations  would  allow  different 
views  to  be  made  of  the  object.  With  an  automated  coded-aperture,  hundreds  of 
different  mask  configurations  could  be  used  in  short  order.  Though  each  configura¬ 
tion  would  not  be  fully  independent  from  the  others,  the  sheer  number  of  different 
perspectives  would  certainly  increase  the  resolution  and  accuracy  of  the  imaging  sys¬ 
tem.  The  programmable  collimator  could  be  fabricated  in  a  converging  or  diverging 
pattern  in  order  to  make  a  coded- aperture  that  zooms  in  or  out  on  a  scene,  a  la 
Rothenbush  [58]. 
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Appendix  A.  MATLAB  Code 


Following  is  a  list  of  MATLAB  code  used  for  this  research.  Note:  Indentation  has 
been  removed  for  brevity.  If  the  user  presses  ctrl+a  and  then  ctrl+i  in  Matlab® 
the  indentation  will  be  automatically  restored. 

1.1  The  Least  Correlated  Mask  Sequence 


Listing  A.l. 

1  7.  Mask  Sequence  that  is  least  correlated  with  the  least  user  ... 
movement  s 

7.  In  this  script,  a  sequence  is  first  constructed  where  each  ... 
element  is  on 

7,  and  off  for  equal  periods  of  time.  The  masks  are  otherwise  ... 
random .  The 

7.  code  will  then  look  at  this  sequence,  and  reorder  the  masks  so  ... 
that  the 

7.  programmable  collimator  user  will  have  the  least  amount  of  ... 
swit  ched 
6  °l  elements 

7.  2d  Lt  Jack  FitzGerald,  AFIT/ENP  ,  180CT11 

clear  all 
close  all 
11  clc 

NumSeq  =  4;  7,number  of  masks  in  the  sequence  MUST  BE  EVEN 
MaskSize  =  10;  7.  length  of  row/column  in  masks  MUST  BE  EVEN,  ... 

SINCE  LOOKING  AT  CENTER  ELEMENT  IN  CORRELATION  MATRIX 
rho.min  =  0.5;  7.  minimum  pinhole  density  (open  elements  divided  by... 
total  elements) 

16  rho.max  =  0.5;  ’/,  max  pinhole  density 

7,7.  Create  a  mask  in  which  all  elements  are  "on"  and  "off"  for  ... 

equal  periods 
7,  of  time 

for  i  =  l:MaskSize 
21  for  j  =  l:MaskSize 

halfcheck  =  0; 
while (half  check  ==  0) 

TimeArray  =  r  ound  ( rand  ( NumSeq  ,  1 ))  ;  7,Make  time  array  to  populate  ... 
sequence 

26  OpenElement  =  0; 

for  time  =  1 : NumSeq 
if (TimeArray (time)  ==  1) 
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OpenElement  =  OpenElement  +  1; 
end 
31  end 

if  ( sum  ( OpenElement )  )  ==  (NumSeq  /  2)  “/.Check  to  see  if  time  ... 

array  is  "on"  and  "off"  for  equal  durations 
half  check  =  1 ; 
end 
36  end 

for  time  =  1 : NumSeq 

Seq  ( t  ime  ,  i  ,  j  )  =  TimeArray  (time)  ;  '/.Populate  sequence  with  time  ... 
arrays 

end 
41  end 
end 

11  Check  correlation  between  one  mask  and  all  others  in  the  ... 
sequence 

46  for  i=l: NumSeq 
for  j  =  1  :  ( i  —  1 ) 

CorrOldMat  =  xcorr2 (reshape (Seq(i , : ,:),MaskSize,MaskSize),reshape( 
Seq(j , : , ,MaskSize ,MaskSize)) ; 

Corr01d(i,j)  =  CorrOldMat (MaskSize , MaskSize ) ; 

CorrOldCoeff (i ,  j  )  =  sum ( sum (corrcoef (reshape (Seq(i,  :  ,  :)  .MaskSize  , . 
MaskSize) ,reshape(Seq(j , : , :) .MaskSize .MaskSize)))) ; 

51  end 

for  j  =  (i+1) : NumSeq 

CorrOldMat  =  xcorr2 (reshape (Seq(i , : , :), MaskSize, MaskSize), reshape ( 
Seq(j , : , :) .MaskSize .MaskSize)) ; 

CorrOld(i.j)  =  CorrOldMat (MaskSize , MaskSize ) ; 

CorrOldCoeff (i ,j)  =  sum ( sum (corrcoef (reshape (Seq(i,  :  ,  :)  .MaskSize  , . 
MaskSize) ,reshape(Seq(j , : , :) .MaskSize .MaskSize)))) ; 

56  end 
end 

OldSeq  =  Seq; 

CorrOldTot  =  sum ( sum ( Corr Old ) )  ; 

61  CorrOldCoeff Tot  =  sum ( sum (CorrOldCoeff ))  ; 

“/.Find  time  open  and  closed  for  average  element 
OpenElements  =  zero s ( MaskS ize ) ; 
for  time  =  1 : NumSeq 
66  for  row  =  1 : MaskSize 
for  col  =  1 : MaskSize 
if (OldSeq (time , row , col) == 1 ) 

OpenElement s ( row , col )  =  OpenElement s ( row , col )  +  1; 

end 
71  end 
end 
end 
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OpenElements  =  OpenElements  ./  NumSeq; 

76  7o70  Reorder  the  sequence  so  that  the  user  has  the  least  number  of  . 
element  changes 

Permut  =  perms  ( 1  :  NumSeq)  ;  '/.Make  matrix  of  all  permutations 
NumPerms  =  numel  ( Permut  (:,  1 ))  ;  7,  Find  number  of  permutations 

81  for  i  =  1 : NumPerms 
changes  =  0; 
for  j  =  l:(NumSeq-l) 
for  row  =  l:MaskSize 

for  col  =  l:MaskSize  7.  Check  to  see  if  each  element  changes  from  . 
one  mask  to  the  next  in  the  seq 

86  if  ((  DldSeq  (  Permut  (  i  ,  j  ),  row  ,  col )  ==  1)  kk  (  DldSeq  (  Permut  (i  ,  j  +1)  . 

row , col ) =  =  0)  ) 
changes  =  changes  +  1; 

elseif  (( OldSeq ( Permut ( i , j ), row , col )  ==  0)  kk  ( OldSeq ( Permut ( i , j .  . 

+  1)  , row , col) ==1)  ) 
changes  =  changes  +  1; 
end 
91  end 
end 
end 

ChangeTr acker ( i )  =  changes; 

end 
96 

[MinChanges  ,  MinChangelnd]  =  min  (  ChangeTracker  )  ;  7.  Find  permutation 
with  least  element  changes 

7.  Create  the  sequence  with  the  least  number  of  element  changes 
for  i  =  1 : NumSeq 

101LeastChangeSeq(i,:,:)  =  reshape(01dSeq( Permut (MinChangelnd , i )  ,  1 :  . . 
MaskSize ,l:MaskSize) .MaskSize .MaskSize) ; 

end 

1.2  ML-EM 


Listing  A. 2. 

7,  Programmable  Collimator  as  a  Coded  Aperture  -  Forward  Problem  . 
Simulation 

2  7.  Use  same  coordinate  system  as  the  MCST  work 
7.  2d  Lt  Jack  FitzGerald,  AFIT/ENP  ,  05DEC11 
7.  Outputs  the  A  parameter  for  the  ML-EM  algorithm 

7.  Updated  on  19DEC11  to  over-resolve  the  detector,  to  account  for 
7  7.  vignetting 

7«  For  use  with  coarser  simulations 
7«  Outputs  A  to  a  text  file 

12  7c  All  distances  in  cm 
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7.  Source  is  in  "front", 


detector  is  "rear 


clear  all 
close  all 
17  clc 

tic 


7,7.  Parameters 

22  7.  DetResFact  =  the  resolution  of  the  simulated  detector  plane,  D, 
relat ive 

7.  to  the  actual  detector  subpixel  resolution 

DetResFact  =  10; 

DetUnitFact  =  20;  7,cm  to  0.5  mm 

7,  Assume  the  detector  is  not  crooked  relative  to  crystal 
27  ColXPos  =  4.54;  7.  Distance  from  front  of  crystal  to  back  of  ... 

collimator 

ColZPos  =  1.5;  7,  Distance  from  bottom  of  crystal  to  bottom  of  ... 
collimator 

ColYPos  =  1.5; 

DetZDim  =  8;  7.  Crystal  is  8  cm  tall 

DetYDim  =  8;  7,  Crystal  is  8  cm  wide 
32  SourceXDist  =  52.6;  7,  distance  of  source  plane  from  detector  ... 

plane 

SourceYDist  =  -1;  7. 

SourceZDist  =  -1;  7.  From  Bottom  of  detector  crystal  to  bottom  of  . 
source  plane 

SourceZDim  =  10;  7.  Source  plane  height 

SourceYDim  =  10;  7.  Source  plane  width 
37  Sour c eRe sFact  =  5;  7.  change  from  cm  to  other  resolution 

InitCounts  =  1000000;  7,Counts  from  source,  isotropic,  arbitrary  . 
period  of  time 

SteelAtten  =  0.93;  7.  Attenuation  through  steel  chamber  dividers 

(100  keV  gamma  through  0.05  cm  of  steel) 


42  7,  Input  the  mask  pattern 
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1 

1 

1 

1] 

G  (10 

:)  = 

[1 

0 

0 

1 

0  0 

0 

0  0 

0] 

7.} 


7.  Save  parameters  as  a  structure 
57  7.  Removed  for  brevity 
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7.7. 

7.  Discretize  the  detector  plane,  D. 

7.  1600x1600  =  ten  times  actual  resolution 

62  7.  DetResFact  =  the  resolution  of  the  simulated  detector  plane,  D, 
relat ive 

7.  to  the  actual  detector  subpixel  resolution 

D  =  zeros (DetZDim*DetUnitFact *DetResFact  , DetYDim*DetUnitFact * .  .  . 
DetResFact) ; 

S  =  zeros ( Sour ceZDim *  Sour ceRe sFact  , Sour ce YD im *  Sour c eRe sFact )  ; 

67 

troubcount  =  0; 

NumSrcRow  =  numel (S ( :  ,  1) )  ; 

NumSrcCol  =  numel (S (1 , : ) ) ; 

72  NumDetRow  =  numel (D  (  :  ,  1) )  ; 

NumDetCol  =  numel (D ( 1  ,:))  ; 

A_NParam  =  25600; 

A_MParam  =  NumSrcRow*NumSrcCol ; 

ForwProbParameters . NumSrcRow  =  numel (S ( : , 1) ) ; 

77  ForwProbParameters . NumSrcCol  =  numel (S ( 1 ,:)) ; 

ForwProbParameters . NumDetRow  =  numel (D ( : , 1) ) ; 

ForwProbParameters . NumDetCol  =  numel (D (1  ,:))  ; 

ForwProbParameters . A.NParam  =  25600; 

ForwProbParameters . A.MParam  =  NumSrcRow*NumSrcCol ; 

82 

clear  D  S 
today  =  date ; 

filename  =  [date  ’ _MLEMForwSim_SrcX ’  num2str ( SourceXDist )  ’ _SrcRes 

’  num2str ( SourceResFact )  ’_Atten’  num2st r ( St  eel At t en *  100 )  ’ _N ’ 

num2str ( A_NParam)  ’ _M ’  num2str ( A_MParam )  ’ MaskKum ’  num2str(... 

MaskNum ) ] ; 

87 

f pr intf ( 1 Comput ing  Forward  Problem  simulation\n\n 1 ) 

A  =  zeros (A_NParam , A_MParam ,  ’ single  1  )  ; 

92  for  SrcCol  =  l:NumSrcCol 
for  SrcRow  =  1: NumSrcRow 

CurrSrcZPos  =  SourceZDist+SourceZDim-(SrcRow -1) /SourceResFact -1/. . 

SourceResFact /2 ;  7»Find  the  current  source  position,  in  cm 

CurrSrcYPos  =  Sour ce YD i st  +  ( Sr cCol - 1 ) / Sour ceResFact + 1 /2/ . .  . 

SourceResFact;  7.  Find  the  current  source  position,  in  cm 

97  DetResponse  =  zeros (NumDetRow , NumDetCol  ,’ single  1 )  ; 

DetResponse2  =  zeros ( 160 , 160 ,  ’ single  1 )  ; 


for  DetCol  =  l:NumDetCol 
102  for  DetRow  =  1: NumDetRow 
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CurrDetZPos  =  DetZDim  -  ( DetRow - 1 ) / DetUnit Fact / DetRe sFact  -  1/2/.. 
DetUni t F  act /DetRes Fact ; 

CurrDetYPos  =  (DetCol-l)/DetUnitFact/DetResFact+l/2/DetUnitFact/. . 
DetRe  sFact ; 

107 

DetToSrcVertAngle  =  atan2 ((CurrSrcZPos -CurrDetZPos) , SourceXDist ) ; 
DetToSrcHor izAngle  =  atan2 ((CurrSrcYPos -CurrDetYPos) .SourceXDist) ; 
CollFrontPlaneVert  =  CurrSrcZPos  -  t an ( DetToSrc Vert Angle )*(..  . 
SourceXDist -ColXPos -5) ; 

CollFrontPlaneHor iz  =  CurrSrcYPos  -  tan (DetToSrcHorizAngle ) * ( . .  . 
SourceXDist -ColXPos -5) ; 

112  CollRearPlaneVert  =  CurrSrcZPos  -  t an ( DetToSrc Vert Angle )*(..  . 
SourceXDist -ColXPos ) ; 

CollRearPlaneHoriz  =  CurrSrcYPos  -  t an ( DetToSrcHor izAngle )*(..  . 
SourceXDist -ColXPos ) ; 

CollFrontRow  =  10-floor  ((  CollFrontPlaneVert -ColZPos  )  *2)  ;  7.  ... 

Collimator  row  that  gamma  goes  through,  front 
CollFrontCol  =  ceil  (  (  CollFrontPlaneHor  iz -ColYPos  )  *2)  ;  /.Collimator 
column  that  gamma  goes  through,  front 
117  CollRearRow  =  10  -  f  loor  ((  CollRearPlaneVert  -  ColZPos  )  *2  )  ;  7.  ... 

Collimator  row  that  gamma  goes  through,  rear 
CollRearCol  =  ceil  ((  CollRearPlaneHoriz -ColYPos  )  *2)  ;  /.Collimator  . 
column  that  gamma  goes  through ,  rear 

SameChamber  =  0; 
inside  =  1 ; 

122 

7.  See  if  gamma  travels  through  the  same  chamber,  from  src 
7.  to  det 

if  (CollFrontRow  ==  CollRearRow)  &&  (CollFrontCol  ==  CollRearCol) 
SameChamber  =  1; 

127  end 

7,Check  to  see  if  gamma  is  within  collimator 
if  (CollFrontRow  >  10) 
inside  =  0; 

132  end 

if  (CollFrontRow  <  1) 

inside  =  0; 

end 

if  (CollFrontCol  >  10) 

137  inside  =  0; 
end 

if  (CollFrontCol  <  1) 

inside  =  0; 

end 

142  if  (CollRearRow  >  10) 
inside  =  0; 
end 
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if  (CollRearRow  <  1) 
inside  =  0; 

147  end 

if  (CollRearCol  >  10) 

inside  =  0; 

end 

if  (CollRearCol  <  1) 

152  inside  =  0; 
end 

if  (SameChamber  ==  1)  &&  (inside  ==  1) 

if  G (CollFrontRow , CollFrontCol )  ==  1 

157  CurrSrcDetDist  =  norm ( [CurrSrcZPos -CurrDetZPos  SourceXDist  ... 

CurrSrcYPos  -  CurrDet  YPo  s  ]  )  ;  °/0Dist  b/t  src  pix  and  Det  pix  ,  cm 

DetResponse(DetRow,DetCol)  =  Ini t Count  s/CurrSrcDetDist~2;  %  1/r  ~2 
attenuation  of  radiation 


end 

else 

162  DetResponse (DetRow , DetCol )  =  0; 
end 

if  (SameChamber  ==  0)  &&  (inside  ==  1) 

RowChanges  =  abs (CollFrontRow  -  CollRearRow); 

167  ColChanges  =  abs (CollFrontCol  -  CollRearCol); 

TotalChanges  =  RowChanges  +  ColChanges; 

LowRow  =  min ( [CollFrontRow  CollRearRow]); 

LowCol  =  min ( [CollFrontCol  CollRearCol]); 

RowsAllOpen  =  1; 

172  ColsAllOpen  =  1; 

for  testrow  =  LowRow : (LowRow+RowChanges ) 
if  G (testrow , CollFrontCol )  ==  0 
RowsAllOpen  =  0; 

177  end 
end 

for  tested  =  LowCol : (LowCol+ColChanges ) 
if  G  (CollFrontRow  ,  tested  )  ==  0 
ColsAllOpen  =  0; 

182  end 
end 

if  (TotalChanges  "=  0)  &&  (RowsAllOpen  ==  1)  &&  (ColsAllOpen  ==  1) 
CurrSrcDetDist  =  norm ( [CurrSrcZPos -CurrDetZPos  SourceXDist  ... 
CurrSrcYPos  -  CurrDet  YPo  s  ]  )  ;  °/0Dist  b/t  src  and  Det,  cm 

187  DetResponse (DetRow , DetCol )  =  ( Ini t Count s / Curr Sr cDet Di s t ~ 2 ) * ( . .  . 

SteelAtten~TotalChanges )  ;  ”/0  l/r~2  attenuation  of  radiation 

else 

DetResponse (DetRow , DetCol )  =  0; 
end 

192  end 
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7„  Turn  the  dead  subpixels  separating  strips  to  0 

if  (mod (DetRow/DetResFact  ,10)  ==  0)  II  (mod (DetCol/DetResFact  ,10) 
==  0) 

197  DetResponse (DetRow , DetCol )  =  0; 
end 

“/oDetRow 

end 

202  “/„  DetCol 
end 

for  tempDetRow  =  1:160 
for  tempDetCol  =  1:160 

207  DetResponse2 (tempDetRow , tempDetCol )  =  sum ( sum (DetResponse ((..  . 

DetResFact  *  tempDetRow -(DetResFact -1) )  : (DetResFact ^tempDetRow)  , ( 
DetResFact  *  tempDetCol -(DetResFact -1) )  : (DetResFact  +  tempDetCol) ) ) 
) ; 

end 

end 

clear  DetResponse 

212 

1  Remove  the  corners  of  the  detector  (set  to  zero) 
for  tempDetRow  =  l:(0.5/8*(8  *DetUnitFact ) ) 
for  tempDetCol  =  l:(1.5/8*(8  *DetUnitFact ) ) 

DetResponse2 (tempDetRow , tempDetCol ) =0 ; 

217  end 

for  tempDetCol  =  6.5/8*(8  *DetUnitFact )  : 8/8* (8  *DetUnitFact ) 

DetResponse2 (tempDetRow , tempDetCol ) =0 ; 

end 

end 

222  for  tempDetRow  =  (0.5/8*(8  * DetUnit Fact ) )  :  (  1  /8* ( 8  *DetUnitFact ) ) 

for  tempDetCol  =  l:(l/8*(8  *DetUnitFact ) ) 

DetResponse2 (tempDetRow , tempDetCol ) =0 ; 
end 

for  tempDetCol  =  7/8*(8  *DetUnitFact )  :8/8*(8  *DetUnitFact ) 

227  DetResponse2 (tempDetRow , tempDetCol ) =0 ; 
end 
end 

for  tempDetRow  =  (l/8*(8  *DetUnitFact ) )  :  ( 1  .  5/8* (8  *DetUnitFact ) ) 

for  tempDetCol  =  l:(0.5/8*(8  *DetUnitFact ) ) 

232  DetResponse2 (tempDetRow , tempDetCol ) =0 ; 
end 

for  tempDetCol  =  7.5/8*(8  *DetUnitFact )  :  8/8* (8  *DetUnitFact ) 
DetResponse2 (tempDetRow , tempDetCol ) =0 ; 
end 
237  end 

for  tempDetRow  =  (7.5/8*(8  *DetUnitFact  )  )  :  (8/8* (8  *DetUnitFact  )  ) 

for  tempDetCol  =  l:(1.5/8*(8  *DetUnitFact  )  ) 

DetResponse2 (tempDetRow , tempDetCol ) =0 ; 
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end 

242  for  tempDetCol  =  6.5/8*(8  *DetUnitFact )  :  8/8* (8  *DetUnitFact ) 
DetResponse2 (tempDetRow , tempDetCol ) =0 ; 
end 
end 

for  tempDetRow  =  (7/8*(8  *DetUnitFact ) )  :  (7 . 5/8* (8  *DetUnitFact ) ) 
247  for  tempDetCol  =  l:(l/8*(8  *DetUnitFact ) ) 

DetResponse2 (tempDetRow , tempDetCol ) =0 ; 
end 

for  tempDetCol  =  7/8*(8  *DetUnitFact )  :8/8*(8  *DetUnitFact ) 
DetResponse2 (tempDetRow , tempDetCol ) =0 ; 

252  end 
end 

for  tempDetRow  =  (6.5/8*(8  * DetUnit Fact ) )  :  (7/8* ( 8  *DetUnitFact ) ) 

for  tempDetCol  =  l:(0.5/8*(8  *DetUnitFact ) ) 

DetResponse2 (tempDetRow , tempDetCol ) =0 ; 

257  end 

for  tempDetCol  =  7.5/8*(8  *DetUnitFact )  : 8/8* (8  *DetUnitFact ) 

DetResponse2 (tempDetRow , tempDetCol ) =0 ; 

end 

end 

262 

AColNum  =  (SrcCol -1) *NumSrcRow+SrcRow ; 

A ( :  , AColNum)  =  reshape (DetResponse2  ,25600 , 1)  ; 

clear  DetResponse2 

267 


end 

clc 

f pr intf ( 1 Comput ing  forward  problem  simulation . \n\nPercent  complete... 
:  7„.2f  ’  ,  (SrcCol  /  NumSr  cCol  *  100  )  ) 

272 

end 

Normalize  =  max (max (A) ) ; 

A  =  A . /Normalize ; 

277 

save (filename , ’ A ’ , ’ForwProbParameters ’ , ’ -v7 .3 ’ ) 
toe 


Listing  A. 3. 

1  ML-EM  algorithm  for  the  Programmable  Collimator  as  a  Coded  ... 
Aperture 

•/.  Use  with  output  from  "  ProgCollCodAp_ForwProbSimForMLEM_Try3  .  m" 
1  Include  "RMC_MLEM",  from  Ben  Kowash ,  in  directory 
°/„  2d  Lt  Jack  FitzGerald,  AFIT/ENP  ,  05DEC11 


clear  all 
close  all 
clc 
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10  Numlter  =  100;  7«  Number  of  iterations  in  the  ML-EM  algorithm 

ForwFilename  =  1 15  -  Dec -20  1 1 . .  . 

_MLEMForwSim_SrcX76_SrcRes5_Atten93_N25600_M10000 1 ; 
load (ForwFilename ) ; 

M  =  ForwProbParameters  .  A_MParam  ;  7«  Number  of  source  pixels 
15  N  =  ForwProbParameters  .  A.NParam  ;  7,  Number  of  detector  pixels 

7.  ensure  that  A  is  normalized 
if  max(max(A))  >  1 
norm  =  max(max(A)); 

20  A  =  A . / norm ; 
end 

7.  Load  the  raw  data 

filename  =  ’ 1 1 12 14 _Conf igl _Co57RodCo60Po int _6hr s  1  ; 

25  extension  =  1  . txt  ’  ; 

Raw  =  load ([ f ilename , ext ens ion] ) ; 

Raw  =  rot90 (Raw  ,  3)  ; 

Raw  =  fliplr(Raw); 

30  7«  create  the  ML-EM  variables 

lam  =  ones  (M ,  1)  ;  7«initial  estimate  of  source  plane,  all  ones 
y  =  reshape  (Raw  ,  numel  (Raw)  ,  1)  ;  7«  Raw  data,  lxN  array 

b  =  1;  7»uncorrelated  background  noise,  per  subpixel 

7«a  =  A’  *  ones  (N  ,  1)  ;  7»  Normalization  parameter  for  A 
35 

clc 

fprintf ( ’ Parameters  loaded.  Beginning  ML-EM  iterative  sequence \n\n ..  . 
’  ) 

F igCount  =  1 ; 

40  for  i  =  1: Numlter 

lam  =  RMC_MLEM ( lam ,  A,  y,  b,  a); 

if  mod ( i , round ( Numlter /25) )  ==  0 

fprintf  (’ \nlterative  sequence  7».0f  percent  complete  ’  ,  i/Numlter  .  .  . 
*100)  ; 

45  end 
end 


Listing  A. 4. 

function  lam.new  =  RMC_MLEM ( lam_old ,  A,  y,  b,  a) 

/  Title:  RMC  ML-EM  routine 

7. 

4  7.  By  :  Ben  Kowash 
7.  Date:  15  Oct  07 

7. 

7.  Description:  This  routine  takes  input  data  and  computes  the 
e  st imat  e 
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7.  lambda_new  based  on  the  properties  of  the  system.  The  method 
used  is  a 

9  7.  maximum - 1 ikl ihood  expectation -maximum  iterative  algorithm. 
lam_new  =  lam_old  .*  (A’  *  (y  ./  (A  *  lam_old  +  b) ) )  ./  a; 

1.3  Coded- Aperture  Decorrelation  Code 

Listing  A. 5. 

7.  Coded  Aperture  Decorrelation  Code 

7.  2d  Lt  Jack  FitzGerald,  AFIT/ENP  ,  21N0V11 

5  /  updated  on  02JAN11 

7.  This  version  includes  the  magnification  of  the  shadow 

clear  variables 
close  all 
10  clc 


7.  Input  the  matched  filter 


G  ( 1  , 

)  = 

[0 

1 

0 

1 

1 

1 

1 

0 

1 

1] 

15 

G  (2  , 

)  = 

[1 

0 

0 

0 

0 

1 

0 

0 

1 

1] 

G  (3  , 

)  = 

[0 

1 

1 

0 

0 

1 

0 

0 

0 

1] 

G  (4  , 

)  = 

[0 

1 

0 

0 

1 

0 

1 

1 

1 

0] 

G  (5  , 

)  = 

[1 

0 

1 

1 

0 

1 

1 

0 

0 

0] 

G  (6  , 

)  = 

[1 

1 

0 

1 

0 

0 

1 

0 

0 

1] 

20 

G  (7  , 

)  = 

[0 

0 

0 

0 

0 

0 

1 

0 

0 

0] 

G  (8  , 

)  = 

[1 

1 

1 

1 

0 

0 

1 

1 

0 

0] 

G  (9  , 

)  = 

[1 

0 

1 

1 

0 

0 

0 

1 

1 

0] 

G  (10 

:)  = 

[1 

1 

0 

1 

1 

1] 

25  7. 7.  Try  using  raw,  subpixel  data 

filename  =  ’ 1 1 1 1 2 1 _Co57_Conf ig3_Po s 1 _ 1 200se c  ’  ; 
extension  =  ’  . txt  ’  ; 

7.  load  and  plot  the  raw  data 
30  RawSub  =  load (  [  f ilename  , ext ens i on ] )  ; 

RawSub  =  rot  90 ( RawSub  ,  3)  ; 

RawSub  =  fliplr (RawSub)  ; 

7oincrease  res  for  raw  image 
35  ResFactorR  =  5; 

MaskDetDist  =  8;  7«dist  from  mask  to  detector,  cm 
MaskSrcDist  =  52.6;  7«dist  from  mask  to  source,  cm 
7»MaskDetDist  =  8;  7odist  from  mask  to  detector,  cm 
40  7«MaskSrcDist  =  22.3;  7ndist  from  mask  to  source,  cm 
M  =  (MaskDetDist+MaskSrcDist ) /MaskSrcDist ; 

for  i  =  1 : 160* ResFactorR 
for  j  =  1 : 160* ResFactorR 
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45  RawSub2(i,j)  =  RawSub ( ce il ( i / Re sFact orR ) , ce il ( j / Re sFact orR ) ) ; 
end 
end 

RawSub  =  RawSub2 ; 
clear  RawSub2 
50 

MaskElPix  =  round  (  1 0*M*  Re  sFact  orR )  ;  ‘/.Number  of  detector  pixels  the 
shadow  of  one  mask  element  will  contain 

for  i  =  1 : 10* MaskElPix 
for  j  =  1 : 10* MaskElPix 

55  GSub(i,j)  =  G ( ceil ( i / MaskElP ix ), ceil ( j /MaskElPix )) ; 

if  GSub  (  i  ,  j  )  ==  0 
GSubBr own ( i , j )  =  -1; 
else 

60  GSubBr own ( i  ,  j )  =  1; 
end 
end 
end 

65  GSubCorrRawSub  =  xcorr2 (GSub , RawSub) ; 

GCorrRWidth  =  numel ( GSubCorrRawSub ( 1 ; 
lobesize  =  (( GCorrRWidth ) -numel ( RawSub ( 1 ,:))) /2 ; 

OHat  =  GSubCorrRawSub ( round ( lobe s ize )  : round ( GCorrRWidth - lobe s ize )  , 
round ( lobesize ) : round ( GCorrRWidth - lobe s ize ) ) ; 

OHat  =  OHat  ./  norm ( norm ( OHat ))  ; 

70 

DecorrAxes  =  1 : numel ( OHat ( 1 ; 

DecorrAxes  =  DecorrAxes  ./  ResFactorR  ./  20; 

GCorrRBrown  =  xcorr2 ( GSubBrown , RawSub ) ; 

75  GCorrRBrnWid  =  numel ( GCorrRBrown ( 1  ; 

lobesizeB  =  (( GCorrRBrnWid) -numel (RawSub (1 ,:))) /2 ; 

OHatBrown  =  GCorrRBrown ( round ( lobes izeB ): round ( GCorrRBrnWid . 

lobesizeB)  ,  round ( lobesizeB )  : round ( GCorrRBrnWid - lobes izeB ) )  ; 
OHatBrown  =  OHatBrown  ./  norm ( norm ( OHatBrown )) ; 

1.4  MCNP  Code 

Listing  A. 6. 

1  °i  MCNP  Input  File  Generator 

7.  2d  Lt  Jack  FitzGerald,  27JULY11,  AFIT/ENP/GNE 
7,  Thesis  Work 

7.  Program  to  make  an  MCNP  simulation  of  the  dynamic,  coded  ... 

aperture 
7,  apparatus  . 

6  7.  Updated  19DEC11 

clear  variables 

close  all 

clc 
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11 

16 

21 

26 

31 

36 

41 

46 

51 

56 


filename  =  ’ MCNPSimlnput ’ ; 
extension  =  ’  . txt  1  ; 

fid  =  f open ( [ f i lename  extension] , ’w’) ; 
f print f ( f id ,[ f i lename  ’\n’]); 

NumParticles  =  100000; 

7.  Position  of  source  and  collimator,  using  "conventional"  scheme 
used  thus 
%  far  in  testing 
7.  Units  in  cm 

coll_x  =  15;  7.  from  detector  plane  to  front  of  collimator 

coll_y  =  1.5; 

coll_z  =  1.5; 

src_x_conv  =  75; 

src_y_conv  =  4; 

src_z_conv  =  4; 

det_thick  =  1;  7»  detector  crystal  thickness,  cm 

det_height  =  8; 
det_width  =  8; 

det_pixheight  =  0.05;  7»  detector  pixel  height 

det_pixwidth  =  0.05;  7.  detector  pixel  width 


7»The  mask  configuration 

f 

hard 

coded 

for 

now 

7.1  means  open 

to  gammas 

J 

0 

means 

closed 

to  gammas 

maskconf ig ( 1 , 

)  =  [0 

1 

0 

1 

i 

1 

1 

0 

1 

1] 

maskconf ig (2 , 

)  =  [1 

0 

0 

0 

0 

1 

0 

0 

1 

1] 

maskconf ig (3  , 

)  =  [0 

1 

1 

0 

0 

1 

0 

0 

0 

1] 

maskconf ig (4 , 

)  =  [0 

1 

0 

0 

1 

0 

1 

1 

1 

0] 

maskconf ig (5 , 

)  =  [1 

0 

1 

1 

0 

1 

1 

0 

0 

0] 

maskconf ig (6  , 

)  =  [1 

1 

0 

1 

0 

0 

1 

0 

0 

1] 

maskconf ig (7  , 

)  =  [0 

0 

0 

0 

0 

0 

1 

0 

0 

0] 

maskconf ig (8 , 

)  =  [1 

1 

1 

1 

0 

0 

1 

1 

0 

0] 

maskconf ig (9  , 

)  =  [1 

0 

1 

1 

0 

0 

0 

1 

1 

0] 

maskconf ig ( 10 

:)  =  [1 

1 

1 

0  1 

1 

1  1 

i  i]; 

maskconfig  =  f lipud (maskconf ig) ; 

7.7.  Convert  src  ,  collimator  and  detector  positions 

det_x  =  - coll_y  ; 

det_y  =  - coll_z  ; 

det_z  =  - coll_x  ; 

src_z  =  src_x_conv  -  coll_x; 

src_y  =  src_z_conv  -  coll_z  ; 

src_x  =  src_y_conv  -  coll_y; 

7.7.  Define  the  cells  ... 


7»Regions  of  apparatus  !=  grid  block 

fprintf (fid  ,  ’  1  2  -9.38  -1  5  10  -11  110  -13  IMP:P  =  l\n’); 
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61  fprintf  (fid  ,  ’2  2  -9.38  -1  5  10  -200  100  -110  IMP:P=l\n’); 

fprintf (f id  ,  ’3  2  -9.38  -1  5  210  -11  100  -110  IMP:P=l\n’); 

fprintf (fid  ,  ’4  2  -9.38  -1  5  10  -11  12  -100  IMP:P=l\n’); 

i  =  l ; 

66 

for  j  =  0:9 
for  k  =  0:9 

if  maskconf ig ( j +1 , k+1)  ==  0 

71  mat  =  3 ;  7.mat  er  i  al  number  for  PEEK 
dens  =  -1.32;  7.g/cm~3 
else 

mat  =1;  7.  material  number  for  air 

dens  =  -0.001293;  7.g/cm~3 

76  end 

if  j  ==  9 
if  k  ==  9 

fprintf  (fid ,  ’"/.i'/.i'/.i  7.i  7,1.6f  -7„i  °/.i  107, i  -l7.i  207.  i  -27.  i  IMP:P  =  l\n’ 

81  , i , j  , k , mat , dens  , i  ,( i  +  1 ),  j  ,( j +1)  , k , (k  +  1 ) )  ; 

else 

fprintf  (fid ,  »7.i7.i7.i  7.i  7.1. 6f  -7.i  7.i  107.i  -I7.i  207.  i  -207.  i  IMP:P  =  l\n 

J 

, i , j , k , mat , dens ,i,(i+l),j,(j+l),k,(k+l))  ; 

end 

86  else 

if  k  ==  9 

fprintf  (fid ,  ’7.i7.i7.i  7.i  7.1. 6f  -7.i  7.i  107.i  -107.  i  207.  i  -27.  i  IMP:P  =  l\n 

, i , j , k , mat , dens ,i,(i+l),j,(j+l),k,(k+l))  ; 

else 

91  fprintf  (fid ,  ’7.i7.i7.i  7.i  7.1. 6f  -7.i  7.i  107.i  -107.  i  207.  i  -207.  i  IMP:P  =  1\ 

n  >  ... 

, i , j , k , mat , dens ,i,(i+l),j,(j+l),k,(k+l))  ; 

end 
end 
end 
96  end 


i=2; 

for  j  =  0:9 
101  for  k  =  0:9 

if  maskconf ig ( j +1 , k+ 1 )  ==  0 

mat  =  2;  '/(material  number  for  CerroBEND 
dens  =  -9.38;  7»g  /  cm  ~  3 
106  else 

mat  =  1;  7.  material  number  for  air 

dens  =  -0.001293;  7.g/cm~3 
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end 


111  if  j  ==  9 
if  k  ==  9 

fprintf  (fid ,  '•/.i’/.i'/.i  7.i  7.1. 6f  -'/.i  %i  107,  i  -l'/.i  207.  i  -27, i  IMP:P=l\n’ 

, i , j , k , mat , dens , i , ( i+1) ,  j ,(j+l) ,k, (k+1) ) ; 
else 

116  fprintf  (fid ,  ’7,i7.i7.i  7.i  7,1. 6f  -7,i  7.i  107,i  -l*/.i  207.  i  -207.  i  IMP:P  =  l\n 

J 

, i , j  , k , mat , dens  , i  ,  ( i  +  1)  ,  j  , ( j  +1)  ,k, (k  +  1))  ; 

end 

else 

if  k  ==  9 

121  fprintf  (fid ,  ’7.i7.i7.i  7,i  7,1. 6f  -7,i  7,i  107,i  -107.  i  207.  i  -27.  i  IMP:P  =  l\n 

, i , j , k , mat , dens ,i,(i+l),j,(j+l),k,(k+l))  ; 

else 

fprintf  (fid ,  '7.i7.i7.i  7,i  7.1. 6f  -7.i  7.i  107.i  -107.  i  207.  i  -207.  i  IMP:P  =  1\ 

n  ’  ... 

, i , j , k , mat , dens ,i,(i+l),j,(j+l),k,(k+l))  ; 

126  end 
end 
end 
end 

131  i =3 ; 

for  j  =  0:9 
for  k  =  0:9 

136  if  maskconf ig ( j + 1 , k+ 1 )  ==  0 

mat  =  2;  '/(material  number  for  CerroBEND 

dens  =  -9.38;  7«g  /  cm  ~  3 

else 

mat  =  3;  7,  material  number  for  PEEK 

141  dens  =  -1.32;  7»g/  cm "  3 

end 


if  j  ==  9 
if  k  ==  9 

146  fprintf  (fid ,  ’7.i7.i7.i  7,i  7.1. 6f  -7.i  7,i  lO'/.i  -17.  i  207.  i  -27.  i  IMP:P  =  l\n’ 

, i , j , k , mat , dens ,i,(i+l),j,(j+l),k,(k+l))  ; 

else 

fprintf  (fid ,  ’7.i7.i7.i  7,i  7.1. 6f  -'/.i  7.i  lO'/.i  -I7.i  207.  i  -207.  i  IMP:P  =  l\n 

, i , j , k , mat , dens ,i,(i+l),j,(j+l),k,(k+l))  ; 

151  end 
else 

if  k  ==  9 
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fprintf  (fid  ,  »7.i7.i7.i  7.i  7,1. 6f  -7„i  7,i  10'/,  i  -107.  i  20%i  -27,  i  IMP:P  =  l\n 
) 

, i , j , k , mat , dens ,i,(i+l),j,(j+l),k,(k+l))  ; 

156  else 

fprintf  (fid ,  >7.i7.i7.i  7.i  7.1. 6f  -7.i  7.i  107.i  -107.  i  20%i  -207.  i  IMP:P  =  1\ 

n  ’  ... 

, i , j , k , mat , dens ,i,(i+l),j,(j+l),k,(k+l))  ; 

end 
end 
161  end 
end 


i=4; 

166  for  j  =  0:9 
for  k  =  0:9 

mat  =  3;  '/material  number  for  PEEK 
dens  =  -1.32;  7«g /  cm  ~ 3 

171 

if  j  ==  9 
if  k  ==  9 

fprintf  (fid ,  ’7.i7.i7.i  7.i  7.1. 6f  -7.i  7.i  107.i  -I7.i  207.  i  -27.  i  IMP:P=l\n’ 

, i , j , k , mat , dens ,i,(i+l),j,(j+l),k,(k+l))  ; 

176  else 

fprintf  (fid ,  »7.i7.i7.i  7.i  7.1. 6f  -7.i  7.i  107.i  -l7.i  207.  i  -207.  i  IMP:P  =  l\n 

J 

, i , j  , k , mat , dens , i , ( i  +  1)  , j  , ( j  +1)  ,k, (k  +  1 ) )  ; 

end 

else 

181  if  k  ==  9 

fprintf  (fid ,  ’7.i7.i7.i  7.i  7.1. 6f  -7.i  7.i  107.i  -107.  i  20%i  -27.  i  IMP:P  =  l\n 


, i , j , k , mat , dens , i , ( i+1) ,j ,(j+l) ,k, (k+1)) ; 
else 

fprintf  (fid ,  ’7.i7.i7.i  7.i  7.1. 6f  -7.i  7.i  107.i  -107.  i  20%i  -207.  i  IMP:P  =  1\ 

n  ’  ... 

186  , i , j , k , mat , dens , i , ( i+1) , j ,(j+l) ,k, (k+1)) ; 
end 
end 
end 
end 

191 

mat  =  4;  7,  mat  number  for  Ge 

dens  =  -5.232;  7.g/cm~3 


7,  Define  the  detector  cells 
196  for  col  =  1:160 
f  or  row  =  1:160 

fprintf  (fid  ,  ’7.05.  Of  0  -6  7  7,i  7.i  7.i  7.i  IMP  :  P  =  1  \n  ’  ,  (  (  col  - 1 )  *  160  +  row 
+  10000)  , -(row-1  +  2000)  ,( row - 1  +  200 1 )  , (col -1  +  1000)  , - ( col - 1  + 1 00 1 ) )  ; 
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end 

end 

201 

"/.  Cell  inside  the  universe  sphere 

"/.fprintf  (fid  ,  "997  1  -0.001293  -1999  8  #(-l  5  10  -11  12  -13)  IMP:P. 
=  1  \ n  ’  )  ;  '/.Air 

“/.fprintf  (fid  ,  ’998  1  -0.001293  -1999  -8  #(-6  7  1000  -1160  -2000  ... 
2160)  IMP  :  P  =  l\n  ’  )  ;  "/.  vacuum 

fprintf (fid 997  0  -1999  8  #(-l  5  10  -11  12  -13)  IMP:P=l\n’); 

206  fprintf (fid 998  0  -1999  -8  #(-6  7  1000  -1160  -2000  2160)  IMP:P=1\ 
n  1 )  ; 

‘/.Cell  outside  the  universe  sphere 
fprintf (fid 999  0  1999  IMP:P=0\n’); 

211  fprintf (fid \n ’)  ; 


7,7,  Block  2,  define  the  surfaces 


216  "/.Planes  for  the  PEEK  or  CerroBEND  limits  in  z 
fprintf (fid 1  pz  0\n’); 
fprintf (fid 2  pz  -l\n’); 
fprintf (fid 3  pz  -4.5\n’); 
f pr int f ( f id ,  "  4  pz  -5.5\n’); 

221  fprintf  (fid  ,  ’ 5  pz  -6\n’); 

fprintf  (fid  ,  ’  6  pz  °/,s\n  ’  ,  num2str  (det_z  +  det_thick/2  )); 

fprintf  (fid  ,  ’  7  pz  °/,s\n  ’  ,  num2str  (det_z  -  det_thick/2)  )  ; 

fprintf  (fid  ,  ’  8  pz  °/.s\n  ’  ,  num2str  (det_z  +  det_thick/2  +  0.5)); 

226  "/.Planes  for  the  outer  edges  of  the  apparatus 
fprintf (fid 10  px  -5\n’); 
fprintf (fid  11  px  10\n’); 
fprintf (fid 12  py  -5\n’); 
fprintf (fid 13  py  10\n’); 

231 

"/,  form  the  planes  that  make  the  grid  block 
for  i  =  100:110 

fprintf  (fid,’"/,  i  py  "/,  .lf\n,,i,((i-100)*0.5)); 
end 

236 

for  i  =  200:210 

fprintf  (fid,""/,  i  px  "/,  .lf\n",i,((i-200)*0.5)); 
end 

241  "/.  form  the  planes  of  the  detector 
for  i  =  1000:1160 

fprintf  (fid  ,  ’"/.  i  px  "/„.2f\n"  ,i  ,  (det_x  +  (i-1000)  *det_pixwidth)  )  ; 
end 

246  for  i  =  2000 : 2160 
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fprintf  (fid,’ "/.  i  py  %  .  2f  \n  ’  ,  i  ,  (  det  _y  +  det  _he  ight  -  (i-2000)  *  .  .  . 
det _p ixhe ight ) ) ; 

end 

'/.Define  the  universe  sphere,  radius  300 
251  fprintf  (fid  ,  1  1999  so  200  $universe  sphere\n;); 

fprintf  (fid,  ’  \  n’ )  ; 


256  "/o  Block  3,  Data  Specs... 


f print f ( f id , ; mode  p\n’); 

"/.Define  the  air,  material  1 

f pr int f ( f id  ,  ; ml  7000  -0.78  $  Nitrogen\n  ’  )  ;  ‘/.Nitrogen  ,  by  atom 
fraction 

261  fprintf  (fid  ,  ’  8000  -0.22  $  OxygenXn’);  % Oxygen  ,  by  atom  . 

fraction 


266 


'/.Define  the  CerroBEND  ,  material 
f pr int f ( f id  ,  ;  m2  83000  -0.500  $ 
fraction 

fprintf  (fid  ,  ’  82000  -0.267 

fraction 

fprintf (fid  ,  ’  50000  -0.133 

fprintf (fid  ,  ’  48000  -0.100 

fraction 


2 ,  from  wiki 
bismuth\n  ’  )  ; 

$  oxygen \n ’ )  ; 

$  tin\n ’ )  ;  % 

$  cadmium\n  ’  ) 


,  11  wood  ’  s  metal 11 

"/.bismuth,  by  atom  ... 

°/o  lead  ,  by  atom  ... 

tin,  by  atom  fraction 
;  "/.  cadmium  ,  by  atom 


"/.Define  the  PEEK 

f  pr  int  f  (  f  id  ,  ’  m3  6000  -0.7889  $  Carbon\nJ  )  ;  '/.carbon  ,  by  atom  .., 

fraction 

271  fprintf  (fid  ,  ’  8000  -0.1661  $  Oxygen\n’);  "/.oxygen,  by  atom 

fraction 

fprintf (fid  ,  ’  1000  -0.045  $  Hydr ogen \n  1  )  ; 


"/.Define  the  Germanium 

f  print  f  (  f  id  ,  ’  m4  32000  -1  $  Germanium\n  ’  )  ;  "/.Germanium 

276 

"/.define  the  source 

"/.fprintf  (fid  ,  1  sdef  pos  %f  °/,f  °/,f  par=2  erg=  .  122\n  ’  ,  src_x  ,  src_y  ,  srcz 
);  "/.  point  isotropic 


"/,  source  emitting  cone  into  collimator 
281  fprintf  (fid  ,  ’  sdef  pos  °/,f  "/.f  "/.f  par=2  erg=  .  122\n  ’  ,  src_x  ,  src_y  ,  src_z 
) ; 

vec_x  =  2.5  -  src_x  ; 
vec_y  =  2.5  -  src_y  ; 
vec_z  =  0  -  src_z  ; 

fprintf  (fid  ,  ’  VEC=“/„f  "/.f  °/„f  DIR  =  dl \n  ’  ,  vec_x  ,  vec_y  ,  ve  c_z  )  ; 

286  cone.rad  =  7 ;  '/.cm  of  cone 

theta  =  atan ( cone.rad/ src_z ) ; 
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f  pr  int  f  (  f  id  ,  ’  SI  1  -1  7,f  l\n’  ,  cos  (theta)  )  ; 

fprintf (fid,’ SP1  0  7.f  7.  f\n’,l-2*pi*(l-cos  (theta)) /4/pi, 2*pi*(l-cos 
(theta) ) /4/pi ) ; 

fprintf (f id  ,  ’ SB1  0.  0.  l.\n') I 

291 

'/.define  the  number  of  particles 
fprintf  (fid,  1  nps  '/,i\n’  ,  NumPart  i  cles  )  ; 


'/.kill  photons  with  energy  less  than  117  keV 
296  fprintf  (fid,  ’CUT :p  j  0.117\n’); 

'/.define  the  PTRAC 

fprintf (fid,  ’ PTRAC  FILTER= . 1 17 ,  . 127 , ERG  EVENT  =  SUR  MAX  =  1000000000\n 
’) ; 

fprintf (fid ,  ’  FILE  =  ASC  WRITE  =  P0S  TYPE  =  P  CELL  =  10001  , \n  ’) 


301 

for  i=10001 : 35599 


if  mod  (  i  ,  10)  ==  0 

fprintf (fid, ’ \n’ ) ; 

306  fprintf (fid , ’  ’); 

end 

fprintf  (fid,  ’  7.  i  ,  ’  ,i)  ; 
end 


311  fprintf (fid ,’ 35600  ’)  ; 
clc 

f printf ([’ Finished  writing  to  file:  ’  filename  extension  ’\n’]); 
316  f close (’ all  ’)  ; 


Listing  A. 7. 

7,  PTRAC  Reader  for  Coded  Aperture  Forward  Problem  S imulat ion 
/  2d  Lt  Jack  FitzGerald,  AFIT/ENP ,  19DEC11 

3 

clear  all 
close  all 
clc 

8  fid  =  f open ( ’ ptrac_6 ’ , ’ r ’ ) ; 

HistList  =  zeros  (35600 , 1 )  ; 

for  i  =  1:2571  /.skip  the  header 
line  =  fgetl(fid); 

13  end 

eventcount  =  0; 

while  ischar(line)  /.check  to  see  if  the  line  has  characters 
18  if  st r cmp  (  1  ine  (7  :  1 1 )  ,  1  3000’) 
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eventcount  =  eventcount  +  1; 

CellNum  =  st r 2double ( line  (46 : 5 1 ) )  ;  %  Find  the  cell  number  for 

the  event 

CellList ( eventcount )  =  CellNum; 
end 
23 

if  st r cmp (line (17:21),’  3000’)  ==  1 
for  i  =  1: eventcount 

Hi st Li st ( CellLi st ( i ) )  =  Hi stLi st ( CellLi st ( i ) )  +  1/ eventcount ; 

weight  the  hist  by  number  of  events  in  history 

end 

28  eventcount  =  0; 
clear  CellList 
end 

line  =  fgetl(fid); 
end 

33  f close ( ’ all ’ ) 

HistList2  =  Hi st List ( 1 000 1  :  35600 )  ; 
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